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SYNOPSIS 


Finite difference and finite element methods have been 
in practice for solving differential equations and their 
convergence analysis has proved that accuracy of these methods is 
of finite order. In problem like numerical weather prediction, 
numerical simulations of turbulent flows and other problems where 
high accuracy is desired for complicated solutions, these methods 
do not give very reliable results. Spectral methods have recently 
emerged as a viable alternative to these methods and have 
performed spectacularly well on many problems (e.g. in fluid 
dynamics where large hydrodynamics codes are now regularly used 
to study turbulence and transition, in numerical weather 
prediction and in ocean dynamics). In this thesis an attempt has 
been made to solve initial boundary value problems using spectral 
methods. 

Spectral methods are based on representing the solution 
as the truncated series of orthogonal polynomial in the spatial 
variable and, in principle, are infinite order accurate. 

Current formulation of spectral methods for solving 
initial boundary value problems employ a spectral discretization 
only in space and rely on finite difference techniques for 
advancing in time. As a result the global accuracy of the method 
is reduced to only finite order unless very small time steps are 
used which is not always practicable. 

In the early 1980s Morchoisne [Rech. Aerosp. 1979-5, 
293-3061 DroDosed a method for solving such systems of equations 
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which was spectral in both space and time. However, even though 
his numerical results were impressive, the method has not 
acquired general acceptance so far. One reason for this is that 
it requires considerably more memory than conventional spectral 
methods. Another possible reason for the neglect of this approach 
is that it lacks a theoretical justification. 

Recently, Dutt [Siam J. Numer . Anal., Vol 27, No 4, 
885-903] proposed a method for solving initial boundary value 
problems which is similar to Morchoisne’s approach in that it 
employs a spectral discretization in both space and time, 
henceforth we shall refer to it as the Galerkin-Col location 
method. The method is set in a Galerkin formulation as it seeks 
an approximate solution which minimizes a weighted sum of the 
residuals in a filtered version of the partial differential 
equations and initial and boundary conditions. 

The solution process, however, effectively amounts to 
collocating the filtered version of the partial differential 
equation and initial and boundary conditions at an overdetermined 
set of points. We show in this thesis that the filtering can be 
dispensed with, and it suffices to collocate the partial 
differential equation and initial and boundary conditions at the 
overdetermined set of points. The solution is then obtained by 
finding a least-squares solution to the overdetermined set of 
equations. It has been proved that the solution thus obtained 
converges to the actual solution at a spectral rate of accuracy 
in both space and time. In practice, the huge, full and 
overdetermined set of equations is solved by iterative techniques 


in which a low order finite difference solver is used as a 
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preconditioner . 

Spectral methods give very accurate approximations to 
hyperbolic problems with smooth solutions. The naive use of 
spectral methods on hyperbolic problems with discontinuous 
solutions, however produces oscillatory numerical results. It has 
been known for some time that these oscillations are in 
themselves not insurmountable but contain sufficient information 
to permit reconstruction of the actual solution. 

The observation that the pointwise convergence of a 
high order polynomial approximation to a discontinuous solution 
is very slow'where as the convergence in a weighted mean is very 
fast suggests to post process the solution obtain by standard 
Collocation or Galerkin methods by a local smoothing in order to 
recover spectral accuaracy. Local smoothing will be carried out 
by a convolution in physical space with a localized function and 
hence by a weighted mean which approximates exceedingly well the 
exact value of the solution. From a mathematical point of view, 
the convergence in the mean can be measured in terms of Sobolev 
norm of negative order. It can be shown that the error between 
the computed and exact solution in a negative Sobolev norm decays 
at rate which depends only on the order of the norm. 

In this thesis Galerkin-col location method is used for 
solving hyperbolic initial boundary value problems and periodic 
initial value problems with nonsmooth data. The whole work is 
presented in the form of five chapters. 

Chapter I is devoted to the introduction of 
Galerkin-Col location method and a brief history of spectral 


methods is given. 
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Chapter II aims at the discussion and implementation of 
the Galerkin-col location method to linear first order hyperbolic 
initial boundary value problems with variable coefficients in one 
space dimension using Chebyshev polynomials as trial functions. 
It is established that the filtering can be dispensed with and it 
suffices to collocate the partial differential equation and 
initial and boundary conditions at the over determined set of 
points . 


In chapter III the system of equations obtained using 
spectral method discussed in chapter II is solved using iterative 
techniques. The preconditioned residual minimization method is 
discussed for scalar and for system case. This chapter also deals 
with the numerical treatment of boundary coditions which is 
essential for effective spectral calculations. Finally, in last 
section of the chapter we show how nonlinear problems can be 
solved spectrally using preconditioning. The efficacy of 
preconditioning is shown by computational results. 

The objective of chapter IV is to discuss the Galerkin 
- Collocation method for solving hyperbolic initial boundary 
value problems using Legendre Polynomials as trial functions. 
Numerical results are given for scalar problems. 

Finally, in chapter V we show that if we filter the 
data and solve periodic initial value problems with nonsmooth 
data using Galerkin Collocation method, then we can recover 
pointwise values with spectral accuracy, provided that the actual 
solution is piecewise smooth. 



CHAPTER - I 


General Introduction 

1 . 1 Introduction 

Many important equations of mathematical physics are 
hyperbolic in nature e.g. Euler equations of gas dynamics. 
Maxwell’s equations, equations of Magneto hydrodynamics , the 
classical and elastic wave equations. There is a constant demand 
for efficient algorithms for solving these equations from 
physicists and engineers. In this thesis a recently proposed 
approach Dutt(1990) is used for solving hyperbolic initial 
boundary value problems numerically which is spectral in 
nature and this method is hereafter called as "Galerkin - 
Collocation method" . 

Before discussing Galerkin - Collocation method we need to 
say a few words about spectral methods. Spectral s method may be 
viewed as an extreme development of the class of discretization 
schemes for differential equations known generically as the 
method of weighted residual (MWR) , (Finlayson and Seriven 1966). 
The key elements of the (MWR) are the trial functions (also 
called approximating functions) and test functions (also known as 
weight functions). The trial functions are used as basis 
functions for a truncated series expansion of the solution. The 
test functions are used to ensure that the differential equation 
is satisfied as closely as possible by the truncated series 
expansion. This is achieved by minimizing the residual i.e. the 
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error in the differential equation produced by using the 
truncated expansion instead of exact solution, with respect to a 
suitable norm. 

The choice of test functions distinguishes between the three 
most commonly used spectral schemes, namely the Galerkin, 
Collocation, and Tau version. In the Galerkin approach , the test 
functions are the same as the trial functions. They are, 
therefore, infinitely smooth functions which individually satisfy 
the boundary conditions. The differential equation is enforced by 
requiring that the integral of the residual times each test 
function be zero. In the Collocation approach the test functions 
are translated Dirac delta functions centered at special, 
so-called collocation points. This approach requires the 
differential equation to be satisfied exactly at the collocation 
points. Spectral Tau methods are similar to Galerkin methods in 
the way that the differential equation is enforced. However none 
of the test functions need satisfy the boundary conditions. 
Hence, a supplementary set of equations is used to apply the 
boundary cond i t i ons . 

The spectral techniques applied in practice involve 
discretization of the spatial variable spectrally and use of a 
finite element or finite difference scheme for marching in time, 
this makes the method over all only finite ordered accurate. In 
early 1980’s MorchoisneC 1979 , 1984) proposed a method for solving 
initial boundary value problems (IBVP) which was spectral in both 
space and time. However, even though his numerical results were 
impressive, the method has not acquired general acceptance so 
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far. One reason for this was that it required considerably more 
memory than conventional spectral method. Another possible reason 
for the neglect of this approach is it lacked a theoretical 
justification. 

Recently Dutt(1990) gave an alternative approach for solving 
initial boundary value problems which is intermediate between 
Galerkin and Collocation methods and we refer to it as the 
Galerkin - Collocation method. Galerkin - Collocation method 
treats discretization in a different fashion and achieves 
spectral accuracy in both space and time. Dutt has provided a 
rigorous frame - work for Galerkin - Collocation method. For 
hyperbolic initial boundary value problems, Dutt has proved that 
the method is stable whenever the IBVP is stable and converges to 
the actual solution at a spectral rate of accuracy in both space 
and time. 

Gal erkin-Col locat ion method involves collocating the partial 
differential equation (PDE) and the initial and boundary 
conditions at an overdetermined set of collocation points. The 
approximate solution is a function, belonging to an appropriate 
finite dimensional space, which minimizes a weighted average of 
the residuals at these points. The finite dimensional space is 
usually the space of some Geggenbaur polynomials and the weights 
arise naturally from the Kreiss-Rauch estimate and the 
Gauss-Lobatto quadrature rule for the polynomials chosen. We 
explain how the method applies to a particularly simple 
hyperbolic initial boundary value problem 

(1.1.1a) u t + u x = 0 -lsxsl.tiO, 
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with boundary condition 

(i. 1.1b) u(-l.t) = g(t) t * 0 , 

and initial condition 

(1.1.1c) u(x,0) = f(x) -1 s x s 1 . 

Suppose Pj(x) denotes a particular kind of Geggenbaur polynomial 

N 

of degree i in x. To find out an approximation u to the solution 

of problem* ( 1 . 1 . 1 ) in the space generated by { P^(x) } Q £ i £ N 

N 

and { P.(t) , . ^ XT . Then u can be written as 

j 0 £ j £ N 

vt N N 

u tx , t) = l E a i i P-(x) P.(t) . 
j=0 i=0 1J 1 J 

N 

¥e want to determine { a^ j^N suc ^ that u approximates the 

solution of (1.1.1) spectrally. To do this we do the following. 

Let { x. } ~ . XT and { t. } ~ . XT denote the Gauss-Lobatto 

1 OiliM J U£J£.N 

quadrature points for the polynomials we have chosen. We will 
collocate the PDE at { (x. ,t.) . . v . 

, l j U£i , j£N 

N N 

Suppose f and g are the unique polynomials of degree N such that 
they interpolate the functions f and g at { x,. ^OiiiN an< * 

{ tj respectively. 



f N (X J ) 

/*S 

•i— 1 

X 

it 

0 

£ 

i i 

N 



■ “V 

0 

£ 

j * 

N 

Define the residuals 






p ij 

, ■ N _ N 
= ( u t + U x 

Mxj.t.) 

0 

£ 

i. j 

£ N 


= u N (x i4 0) - 

f N (x 1 ) 

0 

£ 

i £ 

N 

X j 

= u N (-l,t.) 

J 

- **<*!> 

0 

£ 

j * 

N 


Ideally we will like to choose { a. . such that all these 

1 J Oil 1 $ JjSN 
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residuals are zero. In that case we will have to solve an 

overdetermined system of linear equations with a huge and full 

matrix. Instead we minimize certain weighted averages of the 

residuals p. a., t . , the weights being determined by the energy 
1 J 1 J 

estimates and the quadrature rule. 

It has been proved that the minimization problem we intend 
to solve is equivalent to obtaining a least-squares solution to 
the system of equations formed by enforcing the PDE and the 
initial and boundary conditions at an overdetermined set of 
collocation points. Because of this we have to resort to 
iterative techniques to obtain the least-squares solution. 

Spectral methods give very highly accurate approximations to 
hyperbolic problems with smooth solutions. The naive use of 
spectral methods on hyperbolic problems with discontinuous 
solutions, however, produces oscillatory numerical results. The 
oscillations arising directly from the discontinuity have a 
Gibbs-like high frequency character. It has been known for some 
time that these oscillations are in themselves not insurmountable 
but contain sufficient information to permit reconstruction of 
the actual solution. This is achieved by a filtering of the 
computed values. 

A detailed examination of the effect of filtering for linear 
systems of hyperbolic equations with periodic boundary conditions 
and discontinuous initial data was made by Majda et . ai . ( 1978) . 
They showed that for problems in one space dimension, it was 
possible to achieve a convergence rate of infinite order by a 
proper filtering of the initial conditions and also by applying a 
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filtering during derivative evaluations. However, in two space 
dimensions this infinite order of accuracy can be obtained only 
in a domain which excludes the region of influence and this 
region spreads linearly with time. Moreover, it is not clear as 
to how to handle problems where there are discontinuities in the 
forcing function. 

As opposed to global smoothing , one can post-process the 
solution obtained by standard Collocation or Galerkin methods by 
a local smoothing in order to recover spectral accuracy. The idea 
is based on the observation that while the pointwise convergence 
of a high order polynomial approximation to a discontinuous 
solution is very slow, the convergence in a weighted mean is very 
fast. Local smoothing will be carried out by a convolution in 
physical space with a localized function and hence by a weighted 
mean which approximates exceedingly well the exact values of the 
solution. 

From a mathematical point of view , the convergence in the 
mean can be measured in terms of a Sobolev norm of negative 
order. It can be shown that the error between the computed and 
exact solution in a negative Sobolev norm decays at a rate which 
depends only on the order of the norm. The idea was originally 
developed by Abarbanel et.al. (1986) , Gottlieb et.al. (1985) and 
Mercier(1981 ) . In their formulation the approximate solution is 
obtained by first solving a system of ordinary differential 
equations, arising from either the Galerkin or Collocation method 
and then post processing is applied to the computed solution of 
this semi-discrete system of equations. We are not aware as to 
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how this procedure would deal with problems in which there are 
discontinuities in the forcing function also , instead of just in 
the initial data. We show in this thesis that for hyperbolic 
problems with periodic boundary conditions, it is possible to 
recover pointwise values with spectral accuracy using Galerkin - 
Collocation method even when there are discontinuities in the 
initial data and forcing function, as long as the actual solution 
is piecewise smooth. 

In this thesis, we restrict ourselves to hyperbolic 
equations. The theory will be developed for systems in one space 
dimension. The case of periodic initial value problem with 
nonsmooth data is considered in the last chapter. The problem to 
be addressed here is the linear first order Hyperbolic Initial 
Boundary Value Problem with variable coefficients in one space 
dimension. 

Consider the hyperbolic system of equations 

(1.1.2a) Lu(x,t)=F(x,t) OixsI.tiO 

where L is the differential operator given by 

Lu = u - Au - Bu 

t X 

u(x,t), F(x,t) are vector valued functions with n components. 
A(x,t), BCx.t) are (nxn) matrix valued functions in which each 
entry as a function of x and t is smooth. 

The boundary conditions are given by 
(1.1.2b) tfu(0,t) = g(t) t * 0 

(1.1.2c) Pu(l,t) = h(t) t * 0 

M and P are lxn and kxn matrix valued functions which specify the 
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1 and k inflow variables in terms of the (n-1) and (n-k) outflow 
variables at the respective boundaries. Each entry in them as a 
function of t is smooth, g and h are vector valued functions of t 
with 1 and k components respectively. 

The initial condition is given by 

(1.1. 2d) u(x , 0) = f ( x ) 0 s x s 1 

F, f, g and h are smooth and f, g and h are compatible at the 

space-time corner. 

Since L is hyperbolic A has real eigenvalues and a complete 
set of real eigenvectors. In fact, A is similar to a diagonal 
matrix and without loss of generality we can assume that A itself 
is diagonal. We further assume that the boundary is 
noncharacteristic i.e. zero is not an eigenvalue of A at the 
boundary. Henceforth we assume that A can be written in the form 


A = 


A 

0 


A 


0 

II 


A, .< -51 , . 
lxl lxl 

A (n-1 )x(n-l ) >6I (n-l)x(n-l) 


where A* and A** are both diagonal and 5 > 0 for x e [0,11, t 0 
Kreiss established sufficient condition for problem (1.1.2) 
to be well posed. His condition is known as Uniform Kreiss 
Condition. Rauch modified it further to obtain a priori energy 
estimate for this type of problems. He proved that the following 
estimate holds. 

If u is a solution of (1.1.2) 

T 1 1 


1/2 

( J J II u(x,t)|| 2 dxdt j + ( J II U ^ X ' T ^J| 2 dx j 
0 0 0 


1/2 
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1/2 

( J || u(0, t)|| 2 dt J + [ J II u(l,t)|| 2 dt j 
o 0 


1/2 


T 1 


<; C e 


-rjt 


[T J || F(x,t)|| 2 dxdt + J || f(x)|| 2 dx 
0 0 0 


X X 

+ J II g(t)|| 2 dt + J J|h(t )|| 2 dt ] 
0 0 


1/2 


for all n > rj o , rj Q large enough and T > 0 

The constant C is independent of T, F, f, g, h and r). 

This is known as the Kreiss-Rauch estimate. Using this the 
Galerkin-Col location method is developed for hyperbolic initial 
boundary value problems. 


1 . 2 Layout of the thesis 

The objective of the present thesis is to implement the 
Galerkin - Collocation method to hyperbolic problems with 
different trial functions and to show that the computed solution 
converges to the actual solution spectrally. 

The work embodied in this thesis is divided into five 
chapters . 

Chapter I is concerned with a brief history of spectral 
methods and introduction to Galerkin - Collocation method. 

Chapter II aims at the discussion and the implementation of 
the Galerkin - Collocation method for solving initial boundary 
value problems. The trial functions are the Chebyshev 
polynomials. We show the filtering can be dispensed with and it 
suffices to collocate the partial differential equation and 
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initial and boundary condition at the over determined set of 
points . 

Chapter III deals with preconditioning for scalar and system 
cases and nonlinear hyperbolic problems. We give the numerical 
treatment of boundary condition which is essential for effective 
spectral calculation. Computational results are presented to show 
the efficacy of the preconditioning at the end of each section. 

In chapter IV we show how the Galerkin - Collocation method 
applies to hyperbolic initial boundary value problems using 
Legendre polynomials as trial functions. Numerical results are 
given for scalar problems. 

Finally in chapter V we show that if we filter the data and 
solve periodic initial value problems with nonsmooth data using 
the Galerkin - Collocation method, then we can recover pointwise 
values with spectral accuracy, provided that the actual solution 


is piecewise smooth. 



CHAPTER II 


Galerkin - Collocation Method 


2 . 1 Introduction 

In this chapter we implement a spectral method for solving 
initial boundary value problems which is in between the Galerkin 
and Collocation Methods. In this method the partial differential 
equation and initial and boundary conditions are collocated at an 
overdetermined set of points and the approximate solution is 
chosen to be the least - squares solution to this system of 
equations. In the second section we have discussed the method and 
theoretical results. 

2 . 2 Discussion of Method and Theoretical Results 

In this chapter we restrict ourselves to the case of one 
space dimension. The method we describe, however, is applicable 
to any number of space dimensions. 

We shall shift our initial time from t = 0 to t =-l as this 
will considerably simplify our presentation. Consider the 
differential operator 

(2.2.1) LusUj.- Au x - Bu. 

Here u is a vector-valued function with k components and A and B 
are k x k matrix-valued functions which are smooth functions of x 
and t. We assume the system (2.2.1) is strictly hyperbolic. We 
consider the initial boundary value problem 
(2.2.2a) L u(x,t) = F(x,t) 


,for -1 3£ x j£ 1 ,-l £ t £ 1 , 
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with boundary conditions 

(2.2.2b) M u(— 1 , t ) = g(t) , for -1 s t s 1 , and 

(2.2.2c) P u(l,t) = h(t) ,for -1 £ t & 1 , 

and initial condition 

(2. 2. 2d) u(x,-l) = f(x) , f or -1 s; x s 1 . 

If there are Z inflow variables at the boundary x = -1 then M is 
a Z x k matrix valued function which prescribes the Z inflow 
variables in terms of the (k-.£) outflow variables. Similarly if 
there are s inflow variables at the boundary x = 1 then P is a 
s x k matrix - valued function. Both M and P are smooth functions 
of t. We assume that the initial and boundary data f, g, h and 
forcing function F are smooth and satisfy the compatibility 
conditions which must hold at the space-time corners for the 
solution u to be smooth. Finally we assume that the above 
initial boundary value problem (IBVP) satisfies the Uniform 
Kreiss condition. If the Uniform Kreiss condition is satisfied 
then the IBVP is well posed, i.e. the solution u depends 
continuously on its data. More precisely, it has been proved that 
the following estimate 

11 1 1 

J J || u(x,t) ||^dx dt + J || u(-l,t) || 2 dt + J || u(l,t) || 2 dt 

-1 -1 -1 -1 

r 1 2 

(2.2.3) + J || u ( x , 1 ) || ax 

-1 

^ C [ J J M l| 2(ix dt + J || f(x) || 2 dx + J || g(t) || 2 dt 

-1 -1 -1 -1 
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+ J II h <t> || 2 dtl 

-1 

holds, for some positive constant C. Here the norm || || denotes 

the Euclidean norm. 

One final remark we make is that an IBVP that is well posed 
is ’structurally stable’, i.e. if we perturb the coefficients of 
the differential operator and boundary operator by a small amount 
then the perturbed problem continues to remain well posed. This 
property is crucial for proving that the approximate solution we 
obtain by our method converges to the actual solution of the 
IBVP. 

The method which we now describe applies to general 
Geggenbauer polynomials but in this chapter we shall describe it 
for Chebyshev polynomials. We recall that the Chebyshev 
polynomials T.(y) = cos(j cos *(y)) are orthogonal with respect 

J 

to the weight function 

«Cy ) = 

in the interval [-1,1]. 

Let S p ’ q be the set of polynomials w p,q (x,t) of the form 

P Q 

(2.2.4) w P,q (x,t) = l £ a. . T.Cx) T.(t) , 

i=0 j=0 1J 1 J 

with scalar coefficients a... Similarly, we shall denote by 

x j 

(S p,q ) k the set of polynomials w p,q of the form (2.2.4) if the 

coefficients a. . are vectors with k components. Henceforth we 
^ J 

shall assume that there exists a constant X such that i £ £ £ X. 
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We now define an interpolation operator I p,q which takes a 
continuous function r(x,t) defined on [-1 , 1 ]x[-l , 1 ] and projects 
it into S p ’ q . Thus 

q p 

(2.2.5) I p,q r(x, t) = £ E b. . T (x) T.(t) * r PjQ (x J t) 

j=0 i=0 1J 1 J 

is the unique polynomial belonging to S p,q which interpolates 
r(x,t) at the (p+1) x (q+1) points {(xP.t^)} . n . n . 

Here the points 

x p = cos (irc/p) , 0 £ i i p, and 

t*j = cos (jrc/q) , 0 £ j £ q, 

J 


are the Gauss-Lobatto-Chebyshev points. 

In much the same way we can define a one-dimensional 
interpolation operator I which takes a continuous function s(y) 
defined on [-1,1] and projects it into the space of polynomials 
of degree £ Z. Thus 

B Z B 

(2.2.6) I s(y ) = V b. T . (y) =s^(y) 

i =0 1 1 

is the unique polynomial of degree ^ Z which interpolates s(y) at 

Z 

the (£+1) points {y i = cos( ijr/£)} i=0 ^ . 

We can now use these interpolation operators to define a filtered 
version of the differential operator 

L u = u. - A u - Bu . 

X 


Let 


X p ’ q = I p ’ q A , and 

l p ‘ q = jP.Q B 


be the polynomial interpolants of 


the k x k matrix-valued 
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functions A and B. We now define the differential operator 
(2.2.7) L P,q u = u. - A P * q u - B P,q u , 

Z X 

which can be regarded as a perturbed version of the original 
differential operator Lu. 

Similarly, we define 

= I q tf , and 

P q = I q P . 

We now replace the original IBVP by a filtered version : 

(2. 2. 8. a) L p,q u(x,t) = F(x,t) , for -1 i x , t £ 1, 
with boundary conditions 

(2 . 2 . 8 . b) rf 1 u ( — 1 , t ) * g(t) , for -1 a; t £ 1 , 

(2.2.8. c) P q u ( 1 , t ) = h(t) , for -1 s t a; 1 , 

and initial conditions 

(2 . 2 . 8 . d) u(x, “1 ) = f ( x ) , for -1 a; x s 1. 

The above IBVP will be well posed if we choose p and q large 
enough. In fact since (2.2.8) can be regarded as a perturbation 
of (2.2.2) the following energy estimate, (Dutt 1990) 

J J || u p,q (x,t) || 2 dx dt + J || u p,q (-l,t) || 2 dt 
-1 -1 -1 

+ J^ll u p ’ q (l,t) || 2 dt + J 1 ||u p ' q (x,l) || 2 dx 
-1 -1 

a; C £ J J || L p,q u p,q (x,t) || 2 dx dt + J || u p,q (x,-l) || 2 dx 
-1 -1 -1 

+ J^U u p,q (-l,t) || 2 dt + J^l P q u p,q (l,t) || 2 dtJ , 

-1 -1 


(2.2.9) 
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holds for p and q large enough , with some constant C. Henceforth 
we shall let C denote a generic constant. 

From C2.2.9) the inequality 

J J ||u p,q (x, t)|| 2 dxdt s: C £ J J ||L p,q u p,q (x, t ) || 2 (aC x )oC t )dxdt 
-1 -1 -1 -1 

+ J ||u p,q (x,-l)|| 2 <o(x)dx + J u p,q (-l , t )|| 2 o(t )dt 

-1 -1 

(2.2.10) + J 1 ||P q u p,q (l,t)|| 2 w(t)dtj , 

-1 


follows immediately since the weight function w i 1 . 

We wish to find an approximate solution u p,q (x J t) € (S p,q )^ 
to the above IBVP. Notice that if u p,q (x,t) e (S p,q ) k then 

L p ,q u P.q (x , t ) e ( s 2p ' 2q ) k , 

ifi u p,q (-l,t) e ( S 2 V 

F> q u p ' q (l ,t) e ( S 2q ) s , and 

u p ‘ q Cx,-l> e < S p ) k 

and this suggests that we should accordingly filter our data. 

Let 

F 2p ' 2q Cx,t) = I 2p ’ 2q F(x,t> , 


g 2 q (x,t) 

= I 2q g(t) 

s 


E 2 q Ct) 

= I 2q h(t) 

$ 

and 

? 2 p (x) 

= I 2p f 

s 



be filtered representations of the data. If we were to substitute 
our approximate solution into the IBVP the residuals 

p p,q (x,t) = L p,q u p,q (x 5 t) - F 2p,2q Cx,t) , 
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cr q (t) 


u p,q (-l,t) - g 2q (t) 

(2.1.11) 

i7 Q (t) 

P q 

u p ’ q (l,t) - h 2q (t) 


r p ( x ) = 

u p . q 

(x,-l ) - T 2p (x) , 

would, in general 

not be zero. 

We 

would like to choose our 


approximate solution u p,q (x,t) so that it makes these residuals 
as small as possible and for this we need to define a functional 
which will measure the size of the residuals. 

Accordingly we define a functional 

H p,q ( V P ,q ) = J 1 J 1 ||L p,q v p,q (x,t)-F 2p,2q (x,t)|| 2 w(x)G>(t) dxdt 
-1 -1 

+ J llA^yP^C-l , t)-g 2q (t)|| 2 fc>(t)dt 

-1 

+ J^ljpV^U ,t)-E 2q (t)|| 2 <o(t)dt 

-1 

(2.2.12) + J || v p,q (x,-l )-7 2p (x) J| 2 o(x)dx, 

-1 

where 

p q . 

v p,q (x,t) = Y, t b. . T. (x) T.Ct) e ( S p,q ) K . 
i =0 j=0 1J 1 J 

We choose as our approximate solution the unique u p,q € (S p,q ) k 
which minimizes a functional H p,q (v p,q ) over all v p,q , where 
H p,q (v p,q ) is essentially equivalent to .ff p,q (v P)q ). Now we 
observe that 

p p ’ q (x,t) = L p ‘ q v p ’ q (x,t) - F 2p ’ 2q (x,t) e ( S 2p * 2q ) k , 

n q Ct) = « 4 v'' 5 (-l 1 t)-? 2 «(t) e (sV, 

rj q (t) = P q v p,q U,t) - E 2q (t) € (S 2q ) s , and 
i p (x) = v p ‘ q (x.-l) - 7 2p CD € (S 2p ) k , 
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and so we can exactly evaluate the integrals in (2.2.12) by using 
the very highly accurate Gauss quadrature rules. In particular, 
for the Gauss-Lobatto-Chebyshev rule we have that if s(y) is a 


polynomial of degree £ 2N -1 then 

.1 

(2.2.13) 


N s(y^) 

J 


I 3<y) uCy)dy = 5 .L — it 


N 

where the points y . are given by 

J 

N 

yj = cos(rcj/N) i 
N 

and the weights c . are given by 

J 


j=0 Cj 


0 2» j s N, 


N 


N 


= 2 


c . = 1 

J 


if j / 0 or N, 
otherwise . 


However, there is a stronger version of this rule which we use 
for our particular case. Suppose r(y) is a polynomial of degree 

£ N. Then the inequality , (Canuto et.al. 1987) 

. ,, 2, N. i 

J O N r (y.) 1 9 

(2.2.14) J r (y)w(y)dy £ s 2 J r (y) <a(y) dy 

-1 j =0 c j 1 

holds . 

We can therefore replace the functional fT p ’ q (v p ’ q ) we are 
trying to minimize by an equivalent functional 

2 2q 2p II L p ' q v P ' q Cx 2p ' t 2q >- F 2p ■ 2q (x 2p , t 2q ) || 2 
PjQ( v P»Q} = ft j* £ J J ^ l J 1 1 


j=0 i=0 


o? p x c : 2q 
i j 


(2.2.15) 


2q || v P ’ q (-l,t 2q ) - g 2q (t 2q ) II 5 

** jio ^5 


c . 


2q || 1* v P ' q (l,t 2q > - K 2q (t? q ) || : 
ft V' J J 


2q .E 

4 J=0 


q 7 
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2p || v p ' q (x 2p ,-n - f 2p <x 2p > || : 

7t ^ H 1 1 H 


2p Jo 


2p 

°i 


In fact, using (2.2.14) we conclude that 


B P,q (v Piq ) <; H^Cv*'" 1 ) s 4 ^*^(v^ ) ' 1 ). 

We choose as our approximate solution u p,q e ( S p,q )^ which 
minimizes H p,q . 

In other words, our solution u p,q is given by a least- 
squares solution to the overdetermined system of equations 


rP,q,„P»q 


p,q,„p.q- 


71 

“ 2p — 2q 

4p q c p c 4 


1/2 


f (L P.<l U P.1 - F )Cx? p ,t 2q >} = 0, 
0 s i s 2p , 0 i j i 2q, 


(2.2.16) 


^Tcf 


1/2 


j« q <t 2q > u p ’ <1 C-l,t 2< >> - g(t 2q )j-0, 


0 S j s 2q, 


1/2 


2q of 


n 

2p o 2p 


l^U 211 ) u p ' q C-l,t 2q ) - h(t 2q )j = 0, 


0 s j i 2q, 


1/2 


|u p * q ( X 2 p ,-l) - f)(x^ P )j = 0, 


0 i i £ 2p. 


Here, we have used the fact that T^p, 2q (x^ p ,t? q ) = F (x? p ,t^ q ) 
etc. and so can work with point values of the original data. We 
may write the system (2.2.16) in the form 

(2.2.17) D p,q U p,q = Z p,q , 

where D p,q is a X x i> matrix, U p,q is a i> column vector formed by 
concatenating the point values ^U p,q (x p ,t q ) >. and Z p,q is 

a X-column vector with 
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X = kx(2p+l )x(2q+l ) + (£+s )x(2q+l ) + kx(2p+l), and 
v - kx(p+l )x(q+l ) . 

We emphasize that U p,q denotes the p column vector defined 
above and u p,q (x,t) the polynomial belonging to whose 

point values are the components of U p#q . We wish to find a least 
squares solution to the problem (2.2.17). Clearly ,U Pl<1 must 
satisfy the linear system of equations 

[(D p ' q ) T D p - q ] U p ' <l - (D p ' <1 ) T Z p ' q . 

In (Dutt 1990) it has been shown that the matrix (D p,q )^(D p,q ) 
has an inverse for p and q large enough. Hence the solution to 
the minimization problem is unique. 

To store the filtered representations of the coefficient 
matrices A P,q ,B P,q etc. would place a prohibitive overhead on 
memory requirements for realistic problems; there is a way of 
getting around this, however. Instead of solving the system 

(2.2.17) we choose U P,q which is the least squares solution to 
the unfiltered system of equations 

|lLu P ’ q (x? P , t j Q ) - F(Xj P ,tj q )| = 0, 

0 s i <; 2p , 0 s j i 2q, 
«u*>) ; p -0(-i. t f> - 

(2.2.18) 0 s j s 2q, 

(p(t^ q ) u P,q (l,t^ q ) - h(t? q ) 

I « . „ J J 

0 i j i 2q, 
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2p of 


u p,q (x 2p ,-l) - f(x 2p )l = 0 , 


0 £ i ^ 2p, 


as our approximate solution. 

Notice that the system (2.2.18) may be written as 


(2.2.19) 


D P,q yP.q = z P.q 


which is the same as (2.2.17) except that the matrix D p,q has 
been replaced by D p>q , where D p,q may be regarded as a perturbed 
version of D p,q . If A(x,t) is a smooth function then we know that 

|A(x 2p ,t 2q ) - A P ’ q (x 2p , t 2q ) | 

is spectrally small for all i and j. Using this we conclude that 
the matrix D p,q differs from D P * q by a spectrally small amount 
and hence the difference between U p,q and U p,q , the least 
-squares solutions of (2.2.17) and (2.2.19) respectively, is 
spectrally small. 

We make the argument we have outlined above rigorous in the 
following lemmas and theorem. 

Lemma 2.2.1 

Let v p,q belong to the space of polynomials S p,q defined in 
(2.2.4), Then 

' 1 i r 11 

(2.2.20) J J (v p,q ) 2 w(x)u(t)dxdt s Cj J J (v p,q ) 2 dxdt 

- -1 -1 J ^ -1 -1 J 



2 — 2 /« 

where C 1 = (pq) , for any a > 2 . In particular, choosing 


a = 4 we obtain Cj = E 4 (pq)' 
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Using Holder's inequality we get 
r l l 9 

J J ^ vP,q -' «(x)«(t) dx dt 
^ -1 -1 

r 1 -1 


J J ( ( vP><1 ) 2 j dx dt J J [ «(x)«(t> J dx dt 




t -i -i 


1/0 


where 1/a + 1/0 = 1. Now 


C ( “ <x> P* * C 


dx 


is finite for 1 < 0 < 2 . 

For a fixed value of 0 the right hand side of the last equation 

becomes a constant which we shall denote by D d . Hence we can 

P 


conclude that 

r .1 „1 


J J ^ «(x)o)Ct) J dx dt 


0 


1/0 


- */'■ 


V -1 -l 

Now if sCy) is a polynomial of degree m then by Nikolskii's 
inequality (Canuto et.al. 1987, p 288) 


r 1 r ' 

J f s(y) J dy 

- -1 ^ 


for 1 £ 0 < a . 
Thus we obtain 


1/a 


1/0-1 /a) 


J (.CyljVl 

1-1 J 


1/0 


J 1 / 1 [ (vP ’ q)2 ] 0tdxdt 

“ -1-1 ^ 


1/a 


£ K 2 (4pq) 2Cl_1/a) 


J J ^ vP,<l ^ 2 J dxdt 


^ - 1-1 


Putting E = K 2 (4D_) 2/ ^ we get the required result. 
<x p 
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Lemma 2.2.2 

There are constants I^and K 2 depending on p and q such that the 
estimate 

(2.2.21) Kj || V p ’ q || 2 S || D p ' q V p ' <l || 2 s K 2 || V p ' q || 2 . 

holds. Here 

R 

, and 

, where C denotes a generic constant. 


K x = C/ p‘ 


K 2 = c P 


q p 


Let v p,q (x,t) = Y y a. .T.(x) T.(t) be the polynomial e (S p,q ) l£ 

Li L l » J i J 

j=0 i=0 


such that v P5q (x? , tS) = {V p,q }. . for i=0,...,p ,j=0,...,q. 

1 J 1 i J 


We have that 


II o P,q v p,q II 2 ■ 555 1 l 

11 j=0 i=0 


2 2q 2p || L P '>» v p ’ q C x 2p ,t 2q )|| 2 
_ K V V 1 J 


2p 2q 
c . c . 
i J 


„ 2q || v p ' q C-l,t 2< ») || 2 „ 2q || P q v p ,q ( 1 , t 2< ])|| 2 

2q 2q A 2q 

j=0 c. H j=0 c. M 

J J 

n 2P II v p,q (x 2p ,-l)|| 2 

^ Jo - 


~w 

c i 


Then by (2.2.14) we obtain 


.1 A 


|| D p,q V p ’ q || 2 £, 4 £ J J ||L p,q v p,q ( x,t )|| 2 <a(x)«( t )dxdt 


-1 -1 


+ J Hi* 41 v p,q (-l ,t)|| 2 w(t) dt + J ||P q v p ‘ q (-l,t)|| 2 c>Ct)dt 
-1 -1 


+ J ||v p>q (x,-l)|| 2 o)(x) dx j 
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Now by the inverse inequality for differentiation CCanuto et.al. 
1987, p 295) if v p ’ q e (S p,q ) k then 

J 1 J 1 II §£V P,q (x,t) I f 2 wC x ><oC t >dxdt 
-1 -1 

£ Cp 4 J J II v p,q Cx,t) || 2 w(x)<i)(t )dxdt , 

-1 -1 


and 

J J II §F vP,<1 ^ x,t ^ |J 2 w(x)o(t)dxdt 
-1 -1 

^ Cq 4 J J || v P s q ( x , t ) || 2 o)(x)(o(t )dxdt , 
-1 -1 


where by C we denote a generic constant. 

Further, since sup || A(x,t)|| £ C ,and 

(x , t )e[-l , 1 lx [-1 , 1 ] 

sup || B(x,t)|| £ C , 

(x,t)e[-l,llx[-l,ll 

we obtain 


.1 „1 


J J || L p,q v p,q (x,t) || 2 w(x)«(t)dxdt 


-1 -1 


.1 A 


* C(p 4 +q 4 ) J J || v p,q (x,t> || Z «(x)w(t>dxdt 
-1 -1 


1 H vP ' q(x i- t j ) ll 


q p 


P.q^P + q,.,2 


j=0 i=0 


using (2.2.14). 

Hence we can conclude that 


J 1 J' 1 1 |L P ’ q v p,q (x,t)|| 2 w( x ){o( t ) dxdt <; C (p4 ~ q 9 4 ) ||V p,q || 2 
-1 -1 


Now 
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J ||v p ’ q (x,-l)|J 2 «(x) dxsf £ ||v p ’ q (x p ,-l)|| 2 
-1 j=0 

- i i i n- p q ^ 

j=0 i=0 


Hence we can conclude that 


|'|| vP^Cx.-OH 2 w<x)dx x SE H V P ’ q || 2 , 


and by similar arguments that 


J'll vP'Oc-i,t>n 2 u( »dt x £a || VP-® || 2 , 


J"ll v p,q (i ,t )|| 2 U (t)dt s a || v P’ q || 2 . 


Combining all the above inequalities we conclude that 

II dp-V^h 2 * c ( filial , || v P' q || 2 . 

Now using the condition that there is a constant X such that 

X q 

we obtain || DP’ q V p ' q || 2 s Cp 2 || vP’“ || 2 . 

Next, we have to bound || D p,q V p,q || 2 from below . Using 
(2.2.14) we have that 


.1 .1 


||D p ’ q V p ’ q || 2 *J J | |L P ’ q v p ' q (x,t)|| 2 <o(x)<o(t) dxdt 
-1 -1 

+ J H/f* 1 v p,q (-l,t)|| 2 w(t)dt + J ||P q v p ’ q (l ,t)|| 2 «(t)dt 
-1 -1 


+J ||v p,q (x,-l)|| 2 <o(x)dx. 
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And so by (2.2.10) we can conclude that 

IIDP-V-V * C J 1 J 1 ||v P,q <*.t>l| 2 dx dt , 

-1 "I 

and this together with lemma 2.2.1 and (2.2.14) gives us 

||D P,q v P.a H 2 „ _C_ |v p, q 2 

P 3 « 

Theorem 2.2.1 

The difference between U p,q and U p,q , the solutions of equations 
(2.2.17) and (2.2.19) respectively, is spectrally small 
Proof : 

We have that 

uP ,q . (cdP.QjT D P.q } -l (D P.-J)T z P,q and 
yP.q = ((JP.qjT JP.qj-l ( SP.q)T z P,q. 

Hence 

|| U p,q - U p,q || 

(2.2.22) £ ||{(D P,q ) T D P>q }" 1 || ||(D P ’ q ) T - (D P,q ) T || || Z P,q || 

+ ||(( D P,q ) T D P ' q } -1 -{( D P,q ) T D P,q }" 1 || ||(D P ,q ) T ||Z P,q || . 

Now we know that ||D p,q - D p,q || = 0(-i-^ ) for any s > 0 

P 

And in lemma 2.1 we have shown that 

Kj || V p '<» || 2 S || d p • q V p ’ q || 2 s K 2 || VP’ 9 || 2 , 
which implies that 

1/ Kj i || {(D P,q ) T D P,q ) _1 || * 1/ K 2 . 
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Hence 

(2.2.23) ||{(D P,q ) T D P,q }“ 1 || !|(D P,q ) T -(D P,q ) T || ||Z P,q || 

* 1/KjOC i s )||Z P>q || 

P 

We now estimate the second term in the R.H.S of (2.2.22). 


Clearly 


|| (D P,q ) || £ || (D P,q ) || + || (D P,q ) - (D P,q ) || 

£ 2 ^ K„ + 0(i— ) £ 2 , for p 3 q large enough, 

p 

Put 

M = (D p,q ) T D p,q , and 
N = {(D P,q ) T D P>q . 

Let 

4M = (S p ' q ) T 5 p ’ q - <D p ’ q ) T D p ' q 

Then N = M + AM. It is easy to show that || AM || = 0 (— ) , for 

P S 

al 1 s > 0 . Hence 


„ ~i . II m ’ 1 II 2 II ™ II 

II N - M 1 || £ 


,-l 


1 - II M *|| || AM || 


£ 2 || M 1 || 2 || AM ||, for p,q 


large enough. 
So we obtain 


ii n ~‘ - m 1 n £ ~ § °<V 

Kj p 

Thus 


(2.2.24) ||{(D P,q ) T D P,q }" 1 - {(D P,q ) T D P,q } *|| || (D P,q ) T || Z P,q || 
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And 




0(A->|| Z P ’ q || 


K 


||Z p,q || <; sup ||F(x,t)|| + sup ||g(t)|| 

(x,t)e[-l,l]x[-l,l] tet-1,1] 

+ sup | |hC t ) 1 1 + sup ||f(x)|| . 

tet-1,1] xe [ — 1 , 1 ] 


Now combining (2.2.22), (2.2.23) and (2.2.24) we get the result 

|| U P,q - U P,q || ^ 0(-^) » for a11 s > 0. 



CHAPTER III 


Iterative Solution Of Spectral Equations 
For Chebyshev Polynomials 


3 . 1 Introduction 

The system of equations 
(3.1.1) L P ’ q U P ’ q = Z P ’ q 

we get using spectral methods discussed in chapter II is huge, 
full, ill - conditioned and overdetermined. To obtain a 
least-squares solution to this problem we resort to iterative 
techniques. In the first section of this chapter we describe the 
numerical method for scalar problems. In the second section we 
have discussed the method for the systems case and a numerical 
treatment of boundary conditions is dealt with. Finally in 
section 3 we have considered the method for nonlinear problem . 

For notational convenience we drop superscript p,q in 
u p,q (x,t), U p,q and Z p,q so that u p,q (x,t) ■ u(x,t), U p,q * U and 
Z p,q a Z. Recall that 

P Q 

U p,q ( x , t ) = £ i a ?i q T.Cx) T.(t) 

i=0 j=0 1J 1 J 

denotes the approximate solution of (3.1.1) and 

uP .q = JuP.qcxP.tO)) 

l 1 J i 5:8;:::;$ 

is the vector whose components are the (p+1) x (q+1) values of 
u p,q evaluated at the Gauss -Chebyshev Lobatto points. Ve can then 
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write (3.1.1) in the equivalent form 
(3.1.2) L sp U = Z . 


3 . 2 Preconditioning for scalar problems 


A 


least-square 


solution V of 
L Sp U = Z , 


the equation 


is obtained by minimizing the residual 
H(V) = ||L SP V - Z|| 2 * , 

and this suggests that we should seek the solution by using 
preconditioned residual minimization. For this we need to have 
an approximate inverse, which we shall denote by (L v ) , to the 
matrix L sp and we typically use a low order finite difference 
solver for (L ap ) The method can then be described as : 

1) Given the current guess compute the residual 

R (n) = z _ L sp y(n) 

2) Obtain an improvement for by computing 

y (n) _ ( L a P j ~ 1 R (n) 


3) Update the current value of U 


(n) 


by putting 


u(n + D = + o) 


,(n) 


n 


(n) Cii) 

where o> is chosen as that value of w at which H(U +<i>V ) 

n 

achieves its minimum. o> n can be computed using the formula 


( 


R 


(n) 


sp , r (n) 


) 


"n - ^ L S P L SP v (nj j • 


where ( , ) denotes the standard inner product . 
We now explain each step in more detail. 
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/ __ N 

Given the (p+1) x Cq+1) values of U v which are the point 


( (x p ,tS) > we compute 

1 J j • >p 

J j-o, . . . ,q 


values of u 1, ; (x,t) at the points 
C n ) 

T , the coefficients in its representation as a Chebyshev 
series 


(n) 


P q 


(n) 


u VI1 '(x,t) = £ £ y- / T.Cx) T . ( t ) . 

i=0 j=0 1J 1 J 

This can be implemented using either a two dimensional Fast 
Chebyshev Transform or alternatively by matrix multiplications. 
As the details of this are well known (Orszag 1980 and Pulliam 
et.al. 1981) we do not go into it any further. Since we need to 
compute the residuals on a grid with (2p+l )x(2q+l ) points we pad 

r v 

the representation of u^ with zeros as 


(3.2.1) 


( 2p 2q f , 

u nJ (x,t) = £ £ 7 - T T.Cx) T . ( t ) , 

_■ r\ • a A J J. J 


i=0 j=0 


(n) _ 


where = 0 for i>p or j>q. We can now calculate the values of 

u^ n ^(x,t) at the (2p+l )x(2q+l ) points i( x? p , t^ q )l by 

l 1 J 

using an inverse transform or matrix multiplications. It is now 


an easy matter tocompute the residuals 



P (n) = 
K 1 J 

K n) - * ^ n) - qq<n) - . 


o (n) = 
ij 

[ «tj q ) u <n) (-l ,t^ q ) - 

(3.2.2) 

-5 n> ■ 

f PCt 2q > u (n) C-l,t 2q ) - h(t 2q >') 

V. J J J J 

and 

- 

((u Cn) (x 2p ,-l> - f(x 2p >) . 

For the 

scalar problem we 

shall denote the matrices A and B by a 
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and b. 

The differentiations involved in computing C3.2.2) can be 
implemented using matrix multiplications or transform techniques. 
What is important to note is that these computations can be 
speeded up immensely using vector ization, as was pointed out by 

Orszag (1980). Henceforth we shall denote the vector of residuals 

f(n) (n) (n) (n)'* . p (n) 

Ip 9 o 9 T] s x J oy R • 


2) We now seek an iterative improvement to the vector U 
which we denote by the vector V'" 


(n) 


’= (v (n) (x. P ,t q )\ 

1 1 J 


(n) 


with (p+1) x (q+1) components. Let W denote the prolongation 


c 

> 

0 v 

onto 

'V 

the 

grid with 

(2p + 1) X (2q + 

1) points 


J 5i8; : : 

as 

• >p 

• ,q 

described in 

(3.2.1). We can 

write this 

as 






(3.2.3) 



w <n) = 

p 



where P denotes the prolongation operator. 

C n ) 

It is natural to seek v as the solution of the system of 
equations 

(3.2.4) L ap V (n) = L fd P V tn) = B (n) 

where L 1 is a finite difference discretization of the IBVP on 

#• 

the finer mesh. Then we have 

v <n) _ p-1 £ L fd j -1 R (n> ^ 

where P~* should be interpreted as the generalized inverse of the 
operator P. We can write this in two steps as 

Compute 

Cn) 


a) 


vr (n) = ( L fd 


R 
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b) 


Compute 

y(ll) _ p _ 1 yCXl) 


We describe these steps further. 

a) In computing it is important to choose the finite 

difference operator so that is easily invertible and stable. A 
first or second order implicit approximate factorization code 
based on central differencing ideally fulfills all these 
objectives (Pullaim et.al. 1981 & Steger et.al. 1985). In fact , 
many of the general purpose simulation codes in use in research 
and industry utilize just this approach and it should be possible 
to modify them to perform spectral calculations. Further, since 
these codes involve the solution of a set of independent 
tridiagonal or block - tridiagonal matrix solvers the solution 
process can be vectorized. We refer the interested reader to 
(Pullaim et.al. 1981 & Steger et.al. 1985) for details. 

We indicate the equations obtained from the finite 
difference discretization of 

(3.2.5) w+ n) - a(x , t ) w Cn) - b(x,t) w Cn) = p Cll) (x,t) 

at interior points of the space time square using implicit 

central differencing. Here denotes the vector with 

(2p+l) x (2q+l) values iw^ n ^ (x? p , t? q ) = . 

\ 1 J 1J 





which is second order accurate. 


This can be written as 


“i w i,j + >i w i+i,j = p i,j ■ { °i "i-i.j+i - 

( 3 . 2 . 6 ) ♦ Pi w i,j-l + *i W I+1. j+1 } 


_ /_ 1 . 1 1 a ij ~ - f_ 1 . 1 1 a ij 

a i ~ Y" Ax Ax+Ax' j 2 * a i ” Ax Ax 1 +Axj 2 



where 
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The above equation uses information from the 6 point stencil 

shown in Figure 3.2.1. Thus to advance from time level t . . to t. 

j l j 

we have to solve a tridiagonal system. To initialize the 
procedure we impose the initial conditions 

(3.2.7) w (n ^ = r Cn) , 0 £ i i 2p . 

We can impose the boundary conditions either implicitly or 
explicitly. Inflow boundary conditions pose no problem. Thus if 
x = -1 is an inflow boundary for (3.2.5) we simply impose the 
boundary condition 

(3.2.8) wC 2p,j = a j n) * 0 £ j £ 2q. 

If it is an outflow boundary, however, we either impose the 
partial differential equation at the boundary implicitly or use 
extrapolation techniques (Dutt 1990). Our computational results 
show that the implicit treatment is preferable and so we shall 


say a few words about it. 



x 2p x 2p-l 

Fig 3.2.2 

We use the four point stencil shown in Figure 3.2.2 to obtain the 
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equation 



b 2p , j 


W_ . . < + w_ , . 

2p,J+l 2p-l , j 


w 


2p-l , j+1 


= or. 


This is 

of 

the form 

f 





(3.2.9) 

®2p W 2p-1, j +P 

2p W 2p,j = a j~{ 

“2p 

w 0 

2p- 

l,j+l + ^2p W 2p,j+lj ' 

where 

a 2p 

1 

2At 

a 0 . b 0 

2p » J _ 2p , j 

2Ax 4 

a 2p 

= 

1 

2At 

a 2p , j _ b 2p , j 
2Ax 4 


P 2p 

~ ~ XKt + 

a 2p , j _ b 2p,j 
2Ax — 3 — * 

P 2p 

= - 

1 + 
HX 

a 2p,j _ b 2p , j 

2ax 4 

Note, that 

with this 

treatment of 

the 

boundary 

condition our 

system 

of 

equations 

remains tridiagonal, 

and 

with this we 


conclude our discussion of (a). 


Computation for Scalar case 
Example 3.2.1 

u t - (x+e(x,t)) u x = F(x,t), 


-1 <: x s 1 , -1 s t s 1 


F(x,t> = e(x, t)xsin(j^xe t )xy^e t , 


with initial conditions 

TC — | 

u(x,-l) = sin(j£ xe ) , 

and boundary condition 

u(-l,t) = sinCj^- e^) and uCl.t) = sinCj^ e b ). 
The results obtained for three different values of eCx) . 


1) e^Cx, t)=0.5xsin(3n(x+t>) , 
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2) £ 2 (x , t )=0 . 5xsin(4rc(x+t ) ) , 

3) eg(x,t)=0.5xsin(5jr(x+t)) , 


are shown below. 

N - Number of collection points 


N 

0.5xsin(3re(x+t>) 

0.5xsin(4jtCx+t)> 

0.5xsin(5rcCx+t)) 

33 

<2> (5 . 03x1 0” 3 ) 

<26> (2 . 69x1 0~ 10 ) 

<2> (6.69xl0~ 3 ) 

<30> (2 . 32x1 0 -10 ) 

<3> (6.96X10" 3 ) 
<43> (3.75X10 -10 ) 


Table 3.2.1 


The number inside the first bracket gives the iteration number at 
which the residual given in the second bracket is obtained. 

Example 3.2.2 

u. + x u =0 

t JL 

with initial condition 

u(x,-l) = f(x) 

Results obtained for three sets of initial data 

(i) f(x) = sinCyg x) 

(ii) f(x) = sin(^ x) 

(iii) f(x) = sin(y|^j x) 


are shown in table below. N - Number of collocation points 


N 

sinC T7J?y x) 

s inC-g-j x) 

sinCj^ x) 

33 


<2> (4. 65x10“ 8 ) 

< 13> (8 . 58x1 0 _1 4 ) 

<2> (8.86x10” 8 ) 

< 12> (1.75X10 -12 ) 




Table 3.2.2 
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Notice that in Example 3.2.1 u is an inflow variable at 
x = -1 and x = 1 , that is why the value of u is prescribed at 
both the boundaries. In Example 3.2.2 u is an outflow variable at 
x = -1 and x = 1, and so the value of u at both the boundaries is 
obtained by enforcing the partial differential equation there. 
Orszag C1980) had advocated a filter in which the top one third 
of the frequency components of the numerical solution are removed 
and which he has refered to as the two-thirds rule, for the 
preconditioning to be effective. In table 3.2.3 and table 3.2.4 
we , give comparative results for the two - thirds filter and the 
one - half filter, advocated by us. 

u, - a u =0 

L X 

u(x,-l) = f(x) 

Number of collocation points is 33, 

a = x 


f(x) 

1/2 filter 

V 

2/3 filter 

sin TCtJ x 

<2> (1.71x10" 6 ) 
<58> (4 . 46x1 0 _1 4 ) 

<2> (5 . 22xl0" 6 ) 
<65>(3.32xl0" 7 ) 

It 

sm^x 

<2> (5.21X10 -6 ) 
<51 ) (8 . 67xl0~ 13 ) 

<2> (1.55X10" 5 ) 
<60> (8 . 07x1 0" 7 ) 


<3> (6.60x10" 6 ) 

<2> (3 . 03xl0" 5 ) 

S1 V 

<50>(2.81xl0" 12 ) 

<50X1 .75X10" 6 ) 


Table 3.2.3 
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a = -x 


f ( X ) 

1/2 filter 

2/3 filter 

sin TM x 

<2> (2.1 lxio" 8 ) 
<9>(5.01xl0 -14 ) 

<2> (2.09X10" 8 ) 
<10>(6.43xl0" 10 ) 

sinj^x 

<2> (4.65xl0“ 8 ) 

< 13> (8 . 98x1 0 _1 4 ) 

<2> (4 . 57xl0~ 8 ) 
<16X1 .74X10" 8 ) 

. n 

<2> (8.86x10" 8 ) 

<2> (1.16X10" 7 ) 

S 1 Hy-gX 

— 12 

<12>(1. 75x10 ) 

< 15X7 . 29x1 0~ 8 ) 


Table 3.2.4 

It can be seen from Table 3.2.3 and Table 3.2.4 that the one-half 
filter performs better than the two-thirds filter. 

3 . 3 Preconditioning for System case 

Consider the following hyperbolic system 

( 3.3.1a) w. - A w - Bw = F -1 £ x, t £ 1 

If X 


where 

w = ( w 15 w 2 ) , 


A = 

r aj j (x , t) 

a 12 (x,t) " 

, B = 

r b^ j (x, t) 

b 12 (x,t) ^ 


. a 21 (x,t) 

a 22 (x J t) , 


^ b 21 (x, t) 

b 22 (x,t) > 


F = (FjCx.t), F 2 Cx,t)) T . 


We prescribe boundary conditions 
(3.3.1b) M w(-l,t) = g(t) 

(3.3.1c) P w(l ,t) = h(t) 

and initial condition 


(3.3. Id) 


w(x , -1 ) 


f ( X ) . 
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Here M is a Z x 2 and P is a o x 2 matrix , where 0 £ £ ,4 £ 2 . 
£=0 means there are no boundary conditions at x = -1 and o = 0 
that there are no boundary conditions at x = 1 . 

We use the central differenced discretization of 


(3.3.2) 


w< n) - A(x,t)w^ - B(x,t)w n = p^ n ^(x, t) 


at interior points of the space-time square, to advance the 

solution from time level t.^. to t . . 

J+l J 

As we have seen in the scalar case the equation uses 
information from a 6 point stencil. Thus to advance from time 
level t. + . to t. we have to solve a block tridiagonal system. To 
initialize the procedure we impose the initial condition 


w Cn > =r Cn) 
l ,2q l 


0 £ i £ 2q. 


It is evident that block-tridiagonal matrix solver constitutes 
the major portion of the numerical computation of the standard 
implicit algorithm. Equation (3.3.2) produces a 2 x 2 block 
structure for the implicit operator. T.H Pulliam and Chaussee 
(1981) have given an algorithm which transforms the coupled 
system of equations into an uncoupled diagonal form that requires 
considerably less computational work. 

We describe this algorithm in brief for the system case. An 
implicit approximate factorization scheme for the system can be 
written as 


(3.3.3) 


(I — A. . )Aw . = R. . 
i J ox j l j 


Where Aw^ = Wj +1 “ Wj . and A^ = kCx^.tp , 
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Jt 

as we can handle the lower order term explicitly. Here ^ denotes 
the centered difference approximation to the differential 
operator ^ . 

The matrix A^ has a set of real eigenvalues and a complete 
set of eigenvectors , hence a similarity transformation can be 
used to diagonalize A. . 

1 J 


(3.3.4) A. . = T. . A- ■ TT 1 . 

1J 1J 1J 1J 

and so we write (3.3.3) as 

(3.3.5) 


(T. .T. 1 - At T. . A. . T -1 |_ ) Aw. = R. . 
ij ij ij lj lj $x j 1J 


The modified form of the above equation is constructed by moving 

5 


T outside the difference operator ^ 
diagonal form of the algorithm 


. This results in the 


T ij ( 1 - At \j k > T u Aw j * s ij • 

The modification has introduced an error, but T.H Pulliam and 
Chaussee have shown that the error introduced by the 
diagonal ization is first order in time. The new implicit operator 
( I - At A i j is still block tridiagonal, but now blocks are 

diagonal in form so that the operator reduces to two independent 
scalar tridiagonal operators. 

Numerical treatment of Boundary conditions 


A correct treatment of the boundary conditions is essential 
for an effective spectral calculation. If incorrect boundary 
conditions are imposed on the numerical scheme the resulting 
errors will propagate into thecomputat ional domain. If these 
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errors propagateand/or grow sufficiently rapidly, they will 
destroy the solution. 

Since the system is hyperbolic, A has real eigenvalues and 
a complete set of eigenvectors. So there exists a matrix T such 
that TAT * is diagonal. Equation (3.3.1a) can be rewritten as 

T w. - T A T -1 T w - T B T _1 T w = T F ,or 
t x 

w t - A w x - B w = F. 

Here 

w = T w ,with 

w = ( w is w 2 ) , 

A = T A T _1 

5 = T. T -1 - A T T" 1 - T B T~ 1 , and 

v X 

F = T F . 

The variables Wj and w 2 are called characteristic variables. 

Assume w^ is an inflow variable and w 2 is an outflow 
variable at x = -1. Then the boundary operator M would be of the 
form 

(3.3.6) Afw(-1 , t) = w^-l.t) - a(t) w 2 (-l,t) 

where a(t) is a function of t. Hence the boundary condition at 
x = -1 could be written as 

M w(-l,t) = g(t) . 

For the outflow variable we impose the partial differential 
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equation at the boundary implicitly. If we were to impose (3.3.6) 
in the form 

(3.3.7) w^-l.t.) - «(t.) w 0 (-l,t.) = g(t.) 

1 J J Z J J 

the difference equation would no longer decouple into a set of 

tridiagonal equations but instead would become block tridiagonal. 

We can get around this problem by an approximate treatment of 

the boundary condition. In (3.3.7) we approximate the unknown 

value of w„(t.) by using either 
* J 

a) extrapolation ,or 

b) an explicit finite difference discretization of the partial 
differential equation. 

It is easy to show that both these techniques, which we 
describe below , are GKSO stable for a uniform mesh . 


a( i ) Zeroth order extrapolation 

Here we simply put w 2 (-l,tj) = w 2^ -1,t j+i^ 

a( ii ) First order extrapolation 



At 


At' 

<- Ax-* 


' X 2p 

X 2p-1 


fig 3.3.1 


j + 1 


t 


j+2 


for 2q-l 2: j £ 1 


We define 

V- llt 2q-l ) = W 2 ( ' 1 ’ t 2q ) ' and 
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w 2 (-l,tj) = w 2 C-1 » t j+2- > + (At'+At) x 


At 


0 s j <; 2q-2 . 

b) Explicit difference scheme 

Here we use an explicit difference scheme to compute 

w.C-l.t.) from the values of w(-l,t. .) and w(x 0 . t. .). Since 
1 J J + 1 ZP“1 , J+l 

the boundary condition Cb) does not give good results we omit to 
describe it in detail. 


Computational Results for the System case. 


Example 3.3.1 

. f -0 . 5+0 . OlxsinCrex) 

4 = [ 0.5 

and F(x,t) = { f ^ (x , t ) , f 2 Cx , t ) } 
If our characteristic variables 
Uj(x,t) is an inflow variable at 
variable at x = -1 . 

Case I 

Our solution is 


0 . 5+0 .05xcos(rc(x+t)) J * B ” 0 

are u^(x,t) and u 2 (x,t) then 
x = 1 and u 2 (x,t) is an inflow 


u(x,t) 


sin(^ x e^) 
cos(2x - 3t) x e 


t J e = 0.1 

and the boundary conditions are of the form 
u 2 (-l,t) - 2 x sin(t) Uj (-1 , t) = g(t) , 

Uj(l,t) - e* u 2 (l,t) = h(t) 
with initial data 
u(x , -1 ) = f(x) . 
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We omit to write the rather involved expressions for F,g,h,f . 
Case II 
We choose 


uCx , t) 


sinCjg x e t ) 
cos(2x - 3t) x e 


0.1 


Then boundary conditions are 
u 2 (-l,t) - UjC-l.t) = g(t) , 

Uj Cl , t) - u 2 (l ,t) = hCt) 
with initial data 
u(x,-l) = f(x). 


Number of collocation points is 33, 




A 

B 



Iterations 

Error 

Iterations 

Error 

Case 

I 

4 

56 

9.248x10“ 3 
5.276x10 1U 

2 

50 

2.300x10“ 3 
5.529x10 1U 

Case 

II 

5 

57 

8. 868x10“ 3 
4.319x10 

2 

51 

2.890X10“ 3 
2.383x10 


Table 3.3.1 


In A and B the value of the inflow variable at the boundary 
is obtained using boundary conditions a(i) and aCii) respectively. 
Example 3.3.2 

f 0.5 0.0 ) 

_ [ 0.0 -0.5 J 

u^Cx.t) = x 3 + x t 3 + 100xt^+cos( t ) 

A Q 

u 2 (x,t) = t 4 + (x + t) + sinC x ) 


Case I 
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Case II u^(x,t) = cos s i n ^ x + t)J 

4 3 

UgCx.t) = t + (x + t) + sinC x ) 

Case III Uj(x,t) = cos (j 0 q s ^ n ^ x + 

u 2 (x,t) = sin( 2 x - t) 

Number of collocation points is 33 , 




A 

B 



Iterations 

Error 

Iterations 

Error 

Case 

I 

10 

25 

-02 

3.387x10 
2.048x10 u ' 

2 

12 

7.348x10 J:r 
2.078x10 

Case 

II 

45 

99 

7.521X10 7® 
1.267x10 11 

2 

50 

-04 

4.135X10 V* 
9.393x10 

Case 

III 

2 

25 

1.274X10"®° 
6.163x10 U1 

2 

42 

—ns 

8.812x10 T® 

5 . 435x1 0 _iJ 


Table 3.3.2 


In A and B the value of the inflow variable at the boundary 
is obtained using boundary conditions (b) and a(ii) 

respectively. From table (3.2.2) we conclude boundary condition 
(b) gives poor results in general . 

3.4. Preconditioning for Nonlinear Problems 

We now consider the nonlinear IBVP 
(3.4.1a) u t - A(u) u x = F(x,t) , -1 £ x £ 1 , -1 £ t £ 1 , 

with boundary conditions 

(3.4.1b) Mi(-l,t) = g(t) , -1 £ t £ 1 , 

(3.4.1c) Pu( l,t) = h(t) , -1 £ t £ 1 , 

and initial condition 


(3.4. Id) 


u(x,-l) = f(x) 


-1 £ X £ 1 , 
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and we assume that the solution u 
approximate solution u p>q c (S p,q ) k 

2q 2p 


H P 5 q ( v P 1 q ) 



is smooth. We chose as our 
which minimizes 


A(v* 


,q )v P. 


Q - F )(x 2p ,t 2q )|| 2 
^ J 


j=0 i =0 
2q 


(3.4.2) 


||a*v P,q (-l,t 2q ) - g(t 2q )|| 2 

j=0 

2q 


+ £ ||(Pv P ’ q ( 1 , t 2q ) - h(t 2q )|| 2 

j=0 

2p 

+ *_ £ ||(v p ’ q (x 2p ,-l) - f(x 2p )|| 2 

i=0 

over all v p,q e (S p,q )k .Clearly, this gives rise to a nonlinear 
least - squares problems , which we may write as 

(3.4.3) L sp (U)U = Z . 


We solve this nonlinear least - squares problem by the 
preconditioned residual minimization method as before. We outline 
the main steps below : 

1) To get an initial guess for the solution we letV^^ be 
the solution obtained on the finer mesh with (2p+l )x(2q+l ) points 


by using a 

first 

or second order finite 

difference 

solver for 

the 

nonl inear 

IBVP 

(3.4.1) 

. Then we obtain 

from V (0:> 

by 

truncating 

the 

highest 

half of its 

frequency 

components 

as 


before . 
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til 

2) Suppose at the n stage of the iteration we have an 
approximate solution corresponding to u^ n ^(x,t) , We can 
now calculate the residuals 


(3.4.4) 


p P,q (x,t) 

= 

- A(u Cn) )u£ n) - F 2p,2q (x,t) , 

<7 q (t) 

= 

u Cn) (-l ,t) - g 2q (t) 

rj q (t) 

= P q 

u Cn) (l,t) - h 2q (t) 

r p (x) 

= u (n) 

(x,-l) - T 2p (x) 


We wish to find a correction v(x,t) to u (x,t), corresponding to 
the function v(x,t) e (S p,q )^ we have the vector V. Then v(x,t) 
should approximately satisfy 

(3.4.5a) v. - A(u n )v - A(u n ) v = p^ n <!x,t) , -1 £ t £ 1 , 
v XX 

-1 £ X £ 1 , 

(3.4.5b) £ £ Jv n J(-l,t) = a (n:> (t) ,-l £ t £ 1 , 

(3.4.5c) ^ — Jv n J( l,t) = n Cn) Ct) ,-l £ t £ 1 , 


(3.4.5d) v(x,-l) = T^ n ^(x) ,-l £ x £ 1 , 

which is obtained by linearising (3.4.1) about u^ n \ Thus v can 
be obtained as the solution of a linear IBVP . Hence we can use 
the preconditioning techniques already described to obtain a 
correction V. Specifically, let W be the solution obtained on the 
finer mesh by a finite difference solver for (3.4.5), then V is 
obtained by truncating the highest half of the frequency 
components of W. We compute the relaxation factor u n so as to 
minimize the residual 
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2q 2p 


H P,q C(o) 


4 1 1 II- 

j=0 i =0 


v, - <o A(u 


Cn))v x " p (n) )CXj P J tj q )|j 2 


2q 


+ 3qlll(" [i]^ u n] ( - 1 ’ t f ) - < ’ <n)(t f > II 2 

• 3n li — li * 


(3.4.6) 


j=0 

2q 


55.1 ll(“ [£]’„,»)< II s 


j=0 

2p 


%Y. H- v<x f'- i) - x<n)(x f>n 2 • 


We then define 


j=0 


yCn+l ) _ y (n) + w y 

n 


We remark that our numerical experiments indicate that it is 
enough to consider v as an approximate solution to the partial 
differential equation 

v. - A(u^ n ^) v = ,-1 £t£l , -1 £ x z 1 

L X 

along with the initial and boundary conditions (3.4.5b - 3.4.5d) 
to obtain convergence of the numerical scheme . 


Computational Results for the Nonlinear case 
Example 3.4.1 

We consider the Burger’s Equation 
(3.4.7a) 
with 


u t + u u x = F(x,t) , 


u(x,t) = - 10 + sin(j^ x e*) 
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Then we have to impose the boundary condition 
(3.4.7b) u ( — 1 , t ) = g(t) , -1 £ t £ 1 , 

and no boundary condition at x = 1 . 

The initial condition is 
(3.4.7c) u(x , — 1 ) =f(x) . 

Example 3.4.2 

We consider the isentropic Euler equation 

t 

with 

2 tt 

u(x,t) = x + 10 + cos( sin(x+t)) 

2 2 

p(x,t) = x + 2 + sin(2x-t) and y = 2.0 

The flow is supersonic and so we have to specify no boundary 
condition at x = -1 and both p and u at x = 1 . 

Number of collocation points is 33. 
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CHAPTER IV 


Galerkin - Collocation Method 
For Legendre Polynomials 


4. 1 Introduction 

Galerkin - Collocation method can be applied to general 

Geggenbauer polynomials. In chapter II and chapter III we have 
applied this method to Chebyshev polynomials. This chapter deal 
with Legendre polynomials and the results obtained are similar to 
the one obtained in chapter II. Ve now outline the contents of 

this chapter. In section two we discuss the basic properties of 

legendre polynomials and in the third section we have applied the 
Galerkin - Collocation method to Legendre polynomials. The 
computational results for scalar problemsare presented in section 
four . 

4.2 Properties of Legendre polynomial 

We present here a collection of essential formulae for 

Legendre polynomials which are needed in the later part of this 
chapter. For proofs, the reader may refer to Szego (1939). 
Legendre polynomials { L^Cx), k=0,l,...,} are the eigen functions 
of the Strum Liouville problem 

-(P(x)U' (x))' + Q(x)U(x) = X uCx)U(x) 

in the interval (-1,1) with P(x) = 1 - x , Q(x) = 0 and 

w(x) = 1. L,(x) is an even function of x if k is even and is an 

It 

odd function of x if k is odd. 
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If L, (x) is normalized so that L, (1) = 1 then for any k 


. [k/2 ] 

L *' X) = Jo C ' n 


Z k r 2k-2Z n k-2Z 
^Z x 


where [k/2] denote the integral part of k/2. 

Legendre polynomials are orthogonal with respect to the weight 
function 


w(x) = 1. 


I ‘h 


(x))* dx = (k + 4 ) 


© 


-i 


-1 


The expansion of u(x) € L (-1,1) in term of L^’s is 


u(x) = E L k (x) 


k=0 


where u k given by 


=(k+l/2) J u(x) L^Cx) dx , 


-1 


are called Legendre coefficients. 

Legendre-Gauss-Lobbato points defined as 
N N 

x^ = -1, x„ = 1 and x. = ( 1 , 2 , . . . ,N-1 ) are the roots of 
UN j 

t 

polynomial I^(x). We now consider the discrete Legendre 

transforms 

N _ 

u(x) = V u, L.(x) . 
k=0 K K 

- - N 

u(x) is the unique polynomial of degree £ N such that u(x.) = 

• J 

N N N 

u(x.) V x. = 0,1,..., N where x. are Legendre-Gauss-Lobbato 

J J J 

points. u k are called discrete polynomial coefficients of u(x). 


u k 


The inverse relationship is given by 

v k j=0 J k J J 
where w*j and y k are quadrature weights and normalizing factor, 

respectively and are given by the formulas 


N 


N 

w . = 
J 


N(N+1 ) (Lj^CXj)) 4 


V j = 0, 1 .... ,N 
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and y = £ (L, (x*f) 2 w*J Vk = 0,1 N . 

K j=0 k J J 

4 . 3 Discussion of Method and Theoretical Results 

We shift our initial time t = 0, to t = -1 . 

Consider a well posed IBVP 

(4.3. la) Lu(x,t) = F ( x , t ) -lixil.-lst^l 

with boundary conditions 

(4.3.1b) Mi(-l,t) = g(t) -1 £ t s 1 

(4.3.1c) Pu(l,t) = h(t) -1 £ t £ 1 

and initial condition 

(4.3. Id) u(x,-l) = f(x) -1 £ x £ 1 

which is structurally stable. The operator L ,M,P and function 
F(x,t), g(t), h(t) and f(x) have the same meaning as described in 
section 2.2. 

Let S p ’ q be the set of polynomial w P,q (x,t) of the form 

(4.3.2) w P,q (x,t) = V l a. . L.(x) L.(t) 

j=0 i =0 1J J 

with scalar coefficients a^ and L i be the i th Legendre 
polynomial. Similarly we shall define by (S p,q ) k the set of 
polynomials w p,q of the form (4.3.2) if the coefficients a.^ are 
vectors of length k. Henceforth we shall assume that there exist 
a constant X such that 



Now we define the interpolation operator i p,q which takes a 
continuous function r(x,t) defined on [ — 1 , 1 ]x[-l , 1 ] and projects 
it into S p,q . Thus 

I p,q r(x,t) = Y, Y, a i i L .(t) = 
j=0 i =0 1J J 


(4.3.3) 


r P,q (x,t) 
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is the unique polynomial belonging to S P,q which interpolates 

r(x,t) at the (p+l)x(q+l) point {(x p ,t q )J . Here the 

1 J 1=0,1, ... ,p 

j«o, i, . . . ,q 

point x p , t q are Legendre-Gauss-Lobbato points. 

In much the same way we can define interpolation operator 1^ 

for one space dimension which takes a continuous function sCy) 

defined in [-1,1] and project it into the space of polynomial of 

degree s: 1 . Thus 

1 1 -1 
(4.3.4) I s(y) = £ b.L.(y) = s J (y) 

i=0 1 1 

is the unique polynomial of degree £ 1 which interpolates s(y) 
at the (i+1) Legendre-Gauss-Lobbato points {y.}._ n . .. 

Using these interpolation operators we define a filtered 
version of the differential operator 

(4.3.5a) L P,q u(x,t) = (x , t )-A p ’ q u x (x , t )-B P ’ q u(x , t )=F(x , t ) 

with boundary conditions 

(4.3.5b) rf 1 u ( — 1 , t ) = g(t) 

(4.3.5c) P q u ( 1 , t ) = h(t) 

and initial condition 

(4 . 3 . 5d) u ( x , — 1 ) = f ( x ) . 

The above IBVP is well posed if we choose p.q large enough. 
Since (4.3.5) can be regarded as a perturbation of (4.3.1) the 
following energy estimate 

J 1 J 1 || u p,q (x,t) || 2 dx dt + f\\ u p,q (-l,t) || 2 dt 
-1 -l -1 

+ J^U u p,q (l,t) || 2 dt + J 1 ||u p,q (x,l) || 2 dx 
-1 -1 

£ c j" J 1 J 1 || L p,q u p,q (x,t) || 2 dx dt + J^l U p -(X,-1) || 2 dx 
L -1 -1 -1 
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(4.3.6) + J || rf 1 u P,q (-l,t) || 2 dt + P q u P,q (l,t) || 2 dtJ, 

-1 -1 

holds for p and q large enough , with some constant C. Henceforth 
we shall let C denote a generic constant . 


From the above the inequality 

.1 _1 


-1 -1 


J J ||u p,q (x,t)|| 2 dxdt s C £ J J | |L P 1 q u p,q (x, t)|| 2 dxdt 

1 -1 


1 1 

+ J ||u p,q (x,-l)|| 2 dx + J Ui^ u p,q (-l ,t)|| 2 dt 


-1 


-1 


(4.3.7) 

f ol lows . 


J 1 ||P Q u p,q (l ,t)|| 2 dtj , 

-1 


Notice that if 

u p,q (x,t) 

G 

(S P 

L P,q 

u p ’ q (x,t) 

G 

( s : 


u p ’ q (-l,t) 

G 

( s 

P q 

u p,q (l , t) 

G 

( s 


u P,q (x,-l) 

G 

( s 


2q } s 

p.k 


and 


and this suggests that we should accordingly filter our data. 
Let 


2p.2q (xt ) 

= I 2p,2q F(x , t ) , 

g 2q ( x , t ) 

= I 2q g(t) 

E 2 q Ct> 

= I 2q h(t) , and 

r 2p (x) 

= I 2p f 


be filtered representations of the data. If we substitute our 
approximate solution into the IBVP the residuals 

pP' q (x,t> = , 
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(4.3.8) 


cr q (t ) 


n q ( t) 


T P (x) 


Ifi u p ’ q (-l,t) - g 2q (t) 
P q u p ’ q (l,t) - h 2q (t) 


u p,q cx.-l) - f 2p (x) 


in general not be zero. We would like to choose an approximate 
solution u p,q (x,t) so that it makes these residuals as small as 

possible and for this we define a functional which will measure 

* 

the size of the residuals. 

Accordingly we define a functional 

.1 „1 


(4.3.9) 


• q Cv p ’ q > = I I * v p ' q (x.t)-’F 2p ' 2q Cx,t)|| 2 dxdt 

-1 -1 

♦ J‘||M q v p ' q <-l,0-g 2q Ct>|| 2 dt ♦ J 1 ||P q v p ' q O ,t)-K 2q (t)|| 2 dt 

-1 

+ f\\ V P ’ q (x , - i )-T 2 p (x) (| 2 dx, 


-1 


where V P,q (x,t) G (S P,q ) k . 

We choose as our approximate solution the unique u p,q e (S p,q ) k 
which minimizes a functional H p,q (v p,q ) over all v P,q , where 
H p,q (v p,q ) is essentially equivalent to jEf p ’ q ( v p ’ q ) . 

Now we observe that 


(x , t ) 

= L P * q 

v p ’ q (x 

,t) 

- F 2p ’ 2q (x,t) € < S 2p ’ 2q > k , 

<j q (t) 

= i 

, p ’ q (-l 

,t) 

- g 2q (t) e <S 2q )' t , 

n q (t) 

= , 

, p ’ q (l, 

t) 

- I 2q Ct> e CS 2q ) s , and 

T P (x) 

= v p,q 

(x,-l ) 

— 

7 2p (x) £ (S 2p ) k , 


and so we can exactly evaluate the integrals in (4.3.9) by using 
the very highly accurate Gauss quadrature rules. In particular, 
for the Legendre-Gauss-Lobatto rule we have that if s(y) is a 
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polynomial of degree s 2N -1 then 


(4.3.10 ) 


r 1 N 

I s(y) dy = l w , sCy") 

_i J=0 J J 


where the points y. are Legendre-Gauss-Lobbato points and the 

J 

N 

weights Wj are corresponding weights. Moreover, the inequality 


(4.3.11) 


4 j." ls(yN)l 2 * rV<y)dy * £ w" ls(y»>| 2 

j=0 J J J-l j=0 J J 


holds if s(y) is a polynomial of degree £ N. 

We can therefore replace the functional ,ff Psq (v p ’ q ) we are 
trying to minimize by an equivalent functional 

H P,q (v P,q) = 2 £ || L P,q v p -q(i? |, -t 2<I )-r 2p , 2,l Cx?P,t 2 q)||2 w 2p w 2q 

j=0 i=0 i J i J i J 


(4.3.12) 


+ l' || tP v P>q (-l,t 2q ) - g 2q (t 2q ) || 2 w 2q 
j=0 J J J 

+ T II pCl v P,q (l,t 2q ) - h 2q (t 2q ) || 2 w 2q 
j =0 J J J 


+ l II v p,q (x 2p ,-l) - ? 2p (x 2p ) || 2 w 2p 
i=0 1 


In fact, using (4.3.11) we conclude that 

I7 P ’ q ( v p ‘ q ) s H p,q (v p,q ) s 9 # p,q (v p>q ). 

We choose an approximate solution u p,q e (S p,q ) k which minimizes 
H p ’ q . 

In other words, our solution u p,q is given by a 1‘east- 
squares solution to the overdetermined system of equations 

[w 2p w 2q ) 1/2 { L p - u p - - FMx 2p ,t 2<1 )| = 0. 

0 £ i £ 2p , 0 £ j £ 2q, 
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H 4 ) {^tf) u P,q (-l,t 2q ) - g(t 2q )j= 0, 

(4.3. 13) Oi ji 2q , 

( w j Q ] 1/2 { p<1(t j q) “ Piq (-i*tf) - h(t 2q )} = 0, 

0 £ j £ 2q, 

[«?) 1/2 {u-0<xf.-l, - fXx 2p )} = 0, 

0 £ i S 2p. 

Here, we have used the fact that p 2 P» 2 Q ( x 2p ,t 2q ) = F (x? p ,t 2q ) 

l J ^ J 

etc. and so can work with point values of the original data. We 
may write the system C4.3.13) in the form 

(4.3.14) D P ’ q U P ’ q = Z P ’ q , 

In chapter II we have shown the filtering can be dispensed 
with and it suffices to collocate the partial differential 
equation and initial and boundary conditions at the over 
determined set of points. This result hold even in the case of 
Legendre polynomials. Thus the system (4.3.14) may be written as 


(4.3.15) 


pP.q y P,q = z P,q 


where U Pjq is the least - square solution to the unfiltered 
system of equations 

„f) 1/2 {l > q (xf> ,tf) - FCxf ,t 2q >} = 0, 


0 s: i ^ 2p , 0 £ j £ 2q, 


(w*»] 1/2 / «(t 2q ) > q (-l,t 2q ) - g(t 2q )j = 0, 


(4.3.16) 


0 £ j £ 2q, 


H q ) 


m 2q > t 2q ) - h<t 2q 


r>} - »■ 


0 s j i 2q, 
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( W i P ) 72 ju P>q <xf\-l) ~ f)(x^ p >| = 0 5 

0 £ i s 2p. 

It is simple to prove the analog of the lemmas 2.2.1 and 2.2.2 
and of theorem 2.2.1 for the Legendre case. 

4 . 4 Numerical Results 

The following computational results were obtained using the 
preconditioned residual minimization method discussed in 
chapter III. 


Computation for scalar case 
Examp 1 e 4.4.1 

u. - a(x,t) u - b(x,t) u = F(x,t) 

1 x -1 s x st 1, -1 a£ t s 1 


where 

1 C x+ 1 ) 

a(x,t) = x + sin(n(x+t)) , b(x,t) = sin — ^ — * 

F(x,t) =-^ ■^C n e t sin(n(x+t) ) cos(C n xe t )+sin sin(C n xe t ) 

with initial conditions 

u(x,-l) = sin(C n xe 1 ) , 
and boundary conditions 

u(-l,t> = -sinCC^*) 

uCl,t) = sin(C n e t ) . 

The results obtained for three different values of C n are shown 



in table 4.1.1. 
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N - Number of collection points 


N 

C = — 
n 16 

r - n 
n “ 82 

o 

3 

11 

CD 3 
4* 

17 

<2> ( 1 . 05x1 0 -4 ) 

<48> (6.87X10 -9 ) 

<1> (5 . 26x1 0 _4 ) 

<45> (9.27X10 -9 ) 

<1> (2.73x10 4 ) 

<44> (6.1 lxl 0“ 9 ) 


Table 4.1.1 

where the number in the frist brackets denotes the interation 
number at which the error given in the second brackets is 


obtained . 









CHAPTER V 


Spectral Methods for Periodic Initial 
Value Problems with Nonsmooth Data 

5 . 1 Introduction 

In this chapter we consider hyperbolic initial value 
problems subject to periodic boundary conditions with nonsmooth 
data. We show that if we filter the data and solve the problem by 
the Galerkin-Col location method , discussed in chapter II, then 
we can recover pointwise values with spectral accuracy , provided 
that the actual solution is piecewise smooth. For this we have to 
perform a local smoothing of the computed solution. 

We now outline the contents of this chapter. In Section 2 we 
define the Sobolev spaces we shall work in and describe the 
energy estimates in negative Sobolev norms which are needed in 
this chapter. In Section 3 we briefly describe the 
Galerkin-Col location method and prove that the error between the 
approximate solution computed by this method and the actual 
solution in a negative Sobolev norm decays at a rate which 
depends only on the order of the norm. In Section 4 we explain 
the filtering procedure proposed by Abarbanel , Gottlieb and 
Tadmor and show how it can be applied to the approximate solution 
we obtain by the Galerkin-Col location method to recover pointwise 
values of the solution with spectral accuracy. Finally in Section 
5 we present computational results for the proposed method. 
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5 . 2 Energy Estimates for Hyperbolic Initial Value Problems 
with Periodic Boundaries 


We consider hyperbolic initial value problems with periodic 
boundary conditions. Here after x denotes the vector 

x = ( x llX2 , .. ,x d ). 

Let Q = (0,2jt) a be the space domain and J = (-1,1) be the 
time interval we are considering. Consider the IVP 


d 

(5.2.1) Lu = u. - V A. u -B u = F for (x.t) e Q x J 

lr IX- 

i = l 1 


u = f for (x,t) € Q x {-1}. 

p 

Here u is a vector valued function with values in R and A i , B 
are matrix valued functions. Moreover A ^ , B are smooth functions 
of x and t and per iod ic in x. with per i od 2n , for all j = l,..,d, 

J 

and f and F are periodic in each space coordinate with the same 
period but are not necessarily smooth. 

Before we proceed to describe our numerical method and prove 
its convergence we need to review some a priori energy estimates 
which have been proved for the solutions of the system (5.2.1). 
The interested reader is referred to Rauch(1972) and Taylor(1981) 
for details. 

Let u and v be vector valued functions of x and t and 
2rc-periodic in each space direction. Then we denote 


(u.v) = J J u * V dxdt * and 


QxJ 


QxJ 


II u || = (f f lul 2 dxdt ) 

i x i 


1/2 


Here lul denotes the euclidean norm of u if u is a vector and I AI 
denotes the induced matrix norm if A is a matrix. Similarily we 
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denote 


11 u «. ox/ ( n s ,D “ 4 u , 2 d^at ) i 

ftxj | « I +/3£s 

where a = (a^ jctg » • • soc^) is a multi-index and 

D« u = D ai ...D ad u. 
x Xl x d 

In the same way we define 

( u . v ) = J u * 

1 A . . f . 4 " 


/2 


v dx , and 


Qxi±l 


ftx {±1 } 


II u || = ( J | ul 2 dx ) 

Oxftx { ±i } J x{±1 } 


1/2 


Let 


II » II 


s ,Qx{±l 


, = ( I 

■* A . 


I ID® ul 2 dx ) 1/2 


Qx{±l) laUs 


where a is a multi-index as above. 

We can now state the a priori energy estimates conveniently 
in terms of the Sobolev norms we have just defined. Let ip be the 
solution of the hyperbolic IVP with periodic boundary conditions 
(5.2.2a) Lip = <p for (x,t) e ftxj , 

(5.2.2b) ip = 6 for (x,t) e Qx{-1} , 

where <J> and 6 are smooth functions and periodic in space. Then 

for all integers s £ 0 there exists a constant C , which depends 

s 

only on the smoothness properties of , B such that the 

estimate 


(5.2.3) || IP || + II V II 

s , Qx J s , Qx { 1 } 


* C B ( II ♦ II + || e II , ) 

' s.ftxj s,Qx{-l} 


holds. Henceforth we shall use C and C as generic constants. 

S 

Next, we need to state a version of (5.2.3) for negative Sobolev 
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norms . 

Let w be a function of x and t which is periodic in space. 
Let H = { <J> : <p is a smooth function of x and t which is periodic 
in x and has compact support in t }. We define 


|| w || = sup 

" s ' fixJ <J> e H 


I (w ,<J>) | 

ftxj 

in 


Then H 


s ,QxJ 

is defined to be the completion of H with respect to 


-s ,QxJ 

the above norm. Similarily we define 


|| w || = sup 

-s ,Qx{-l } 0 e H 


I (w ,$>) I 

Qxc-n 

Tmi 


s ,ftx{-l } 

With these definitions we can now state the energy estimates in 

"negative" Sobolev norms. For any s £ 0 there exists a constant 

C , which depends only on the smoothness properties of A. , B 
S 1 

such that 


iifii ♦ ii *ii 

-s.ftxJ - s , Qx { 1 } 


(5.2.4) 


S c 3 (II * II + II e II ) , 

^ -s,QxJ -s , Qx {-1 } 


where tp is the solution of (5.2.2), for all 4> and 0. For the sake 
of completeness we shall provide the proof of (5.2.4) below , 
which is very similar to an analogous result proved by 
Rauch( 1972) . 

We consider the following hyperbolic IVP with periodic 
boundary cond i t i ons 

d 

(5.2.5a) L* w = -w t + £ (aT w) x - B t w = x for (x»t) e QxJ , 

i = l i 

(5.2.5b) w = fJ * or € > 


which is the adjoint of (5.2.2). Notice that for this problem we 
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let time run backwards. The following energy estimate is then 
valid for the solution w of the adjoint problem . 

For every s £ 0 there exists a constant C g which depends 
only on the smoothness properties of A. , B such that 


(5.2.6) 


II « II + II w || 

s.QxJ s,Qx{-l} 

* C ( II x II 


s.QxJ 


II HI , ) 

s ,Qx{ 1} 


holds. 


Let ip be the solution of (5.2.2). An integration by parts 
yields 

$ 

(5.2.7) (tp.L w) = (Ltp.w) + (ip, w) - (ip,w) , 

QxJ QxJ Qx{-1} Qx{l} 

since the integrands are periodic in space. 

Let w be the solution of the adjoint IVP with periodic 
boundary conditions 

$ 

(5.2.8a) L w = x for (x,t) e QxJ , 

(5.2.8b) w = 0 for (x,t) e Qx{l} . 

Then by (5.2.7) we have 

(5.2.9) I (ip,x) I * || MP || x || w || 

QxJ -s.QxJ s.QxJ 

+ II V II , x II w II , i • 

-s,ftx{-l} s,Qx{-l} 

But using the estimate (5.2.6) we have 

II W II + II * II * C s II x II . 

s.QxJ s,Qx{-l) s.QxJ 

and this together with (5.2.9) gives 

(5.2.10) l(ip,X> I * C ( || Lip || 

QxJ -s.QxJ 

+ ll V II ' . ) x II X || 

-s,Qx{-l} s.QxJ 

Thus from (5.2.10) we obtain 

<5.2.10 || * II * C ( || L» || ♦ || V || 

-s , fix J -s.QxJ -s ,Qx{-l } 


)• 
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Next , let w' be the solution of the adjoint IVP with periodic 
boundary conditions 

(5.2.12a) Lw=0 for (x,t) e QxJ , 

(5.2.12b) w = n for (x,t) e Qx{l} . 

Then (5.2.7) takes the form 

(5.2.13) (y,/i) = ( Ltp , w ) + (tp,w) 

Qxiu fixj dx{-n 

and by (5.2.6) the estimate 

(5.2.14) 


II w II + II w II - ~ s 

s.QxJ s,Qx{-l} s s j Qx { 1 ) 


* C R II M II 


is valid. 

Now (5.2.13) and (5.2.14) give 




Qx { 1 } a - s , Qx J 


I * C s ( || Ly, || 

+ ll V> || 


) x II HI 

s.Qxl-l} s , Qx { 1 } 


from which we obtain 
(5.2.15) || ip || 

-s,Qx{l} “ -s,QxJ 

Combining (5.2.11) and (5.2.15) we get (5.2.4). 


^ c s ( II L V II 


II II 


-s ,Qx{— 1 } 


) 


5 . 3 Error Estimates for Blended Four ier-Legendre Methods 
for Periodic Problems with Nonsmooth Data 


Henceforth we shall take Q =(0,2it) since the results we 

state carry over tb the general case Q =(0,2n) d in a 

straightforward manner. We now introduce some notation. For each 

N 

integer N we denote by n the space of algebraic polynomials in 

the variable t of degree upto N. For each integer M we denote by 
M 

S the space 

S M = span { e lkx I -M £ k ^ M } . 
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Then we def ine the space V 


M,N 


as the tensor product 


V«' N = { * : *<*,» - l E a mn a 1 ” V » } , 

n=0 m=-M mn n 

where L n (t) is the Legendre polynomial of degree n. Henceforth we 
shall assume that there exists a constant X such that 
1/X i M/ N £ X. 

o 

For any function w periodic in x , which also belongs to L (QxJ) , 
let w denote the projection of w into (V** ,N ) P , i.e. 


d M,N 
P w 


N 

= E 


M 

E 


n=0 m=-M 


xmx j 

w mri e L ( t ) 
mn n 


where 


w 


oo 

E 

n=0 


00 

E 

m=- 


w 


imx 


•oo 


mn 


L n C« 


Henceforth we shall denote P M,N w by w 

The following results are well known , (Canuto et.al. 1988). 
If W € H k,QxJ then 

(5.3.1) || w - w M,N || £ C N _k || w || 

O.QxJ k.QxJ 


Moreover 

(5.3.2) 



II * II w II ■ 

O.QxJ O.QxJ 


Also we have 

(5.3.3) || w - w M,N || * 

1 ,QxJ • 


C N 21 ~ k || w || 

k.QxJ 


for all 0 <. 1 £ k . 

Next , we introduce the norm 

|| w || = max ( ess sup I D“ w I ) . 

s jOO.QxJ a+/3£s (x,t) e QxJ 


Then we have 
(5.3.4) 


ii " M.N ,, 
|| w - w || 


1 jOd.QxJ 


* C N 21 ~ k || w || 


k jOd.QxJ 


for all 0 < 1 s k. 
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If s(x) is a periodic f unct ion belonging to L 2 (Q) we define 

M 


P M ’°s = l s _ e imx 


m=-M 


m 


= s 


M 


oo 


where s(x) = Y s e 
u m 
m=-oo 


imx 


Similarily , if h(t) e L 2 (J) we define 


,0,N 


N 


h = £ h n L n (t) = H * 
n=0 n n 


00 


where h(t)= Y h L(t) 

Lj - n n 


n=0 


We have results similar to (5 . 3 . 1 )-(5 . 3 . 3) for the above. 
Let 

l M-l.N-l = p M-l , N-l A 
| M-l , N-l _ p M-l , N-l B 
p 2M-1 , 2N-1 _ p 2M-l,2N-l p 
j 2M-1 = p 2M-l , 0 f 

We define the differential operator 

L M * N w = w t - A w - B w. 

L X 

We choose as our approximate solution 


M ,N 


M XT N M 

e (v M > N )P = { * : <p(x,t) = l l a 

rs 


n=0 m=-M 


imx T 

> e L ( t ) j 
mn n 


a e R P } 
mm 


which minimizes 


(5.3.5) 38^’ N (w M,N ) = J J IL M,N w M * N - f 2 * 4-1 ’ 2 **" 1 | 2 dxdt 

QxJ 

I w M,N (x t -l) - f 2M-1 (x) I 2 dx 


I 


Qx{- 1 } 

ii M,N ^ ...M.N.p 
over all w € (V ; . 

The above problem reduces to obtaining a 1 east-squares solution to 
an overdetermined set of equations obtained by collocating the 
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modified equation 

l M,N w M,N _ p2M-l , 2N-1 


and the initial conditions at an overdetermined set of points. We 
briefly explain this. 

Let Xj = ni/ M , 0 s i * 2M-1 ,and let lx** } . n M be the 
Gauss-Lobatto-Legendre points with t q = -1 and t n = 1. Notice that 
C L M,N w M,N - f 2M "l»2N-l ) £ ( y 2M-l , 2N-1 } p 


and 

< w M,N (x,-l) - f 2M_1 ) e C S 2M_1 ) P . 
Hence we have that 

2N 4M-1 


S^Cw^N) =1 l a M * N !(L M * N w M ‘ N - F^-^^-hcx^.r 2 ^ 
j=0 i =0 1J 1 J 


(5.3.6) 


•T 1 e" i w M -V M .-n - f 2 M - 1 ( X f)i 2 

i=0 1 1 


M N M 

where a^j and are appropriate constants obtained from the 

Gauss-Lobatto integration formulae. Thus obtaining a solution to 

(5.3.5) is equivalent to solving a least-squares problem. It has 

been shown in chapter II that if we choose our approximate 
_M , N 


solution v 


such that it minimizes the modified functional 


i = 2 £ “- 1 M N |(l „M,N_ p2M-X.2N-l )( 2M 2N, ,2 

j=0 i =0 1J 1 J 


4M-1 


♦ r rf I w M ’ N (xf,-n - f 2 «->cx 2M )i 

i=0 1 


then we would be committing , in addition , only a spectrally 
small further error. There is therefore no need to filter the 
coeffecients A and B in practice. We are interested in another 
aspect of this minimization procedure. Our approximate solution 
v M,N is the unique polynomial belonging to (v M,N ) p which 


satisfies 
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(5.3.7) J J <L M ’ N v M ’ N - p2M-l , 2N-1 } * (L M.N y M, N) dxdt 

QxJ 

J (v M ’ N - f 2M -‘ )* y M ’ N d* = 0 , 

Qx(-l) 

_ M,N ,, 7 M.N.p 

for all y ' e (V ) v . 

We shall now use the above relation to prove that 


U l M,N y M,N _ j;2M-l , 2N-1 || 


-s.QxJ 


£ C N 
s 


1 — s 


and 


II v M ’ N - f 211 - 1 || 


-S ,Qx{-l } 


£ C N 1 S , 
s 


for any s > 1. 

In addition to this we shall also prove 


and 


|| l M,N u _ jp2M— 1 , 2N-1 


II « " ? 2M-1 II 


-s , QxJ 


£ C a N~ S . 
s 


<; C N~ s . 


-s ,Qx {-1 } 

With these results established we can prove Theorem 5.3.1 and the 
reader is advised to proceed directly to the theorem on page — 77 
and continue his perusal of how these result are established only 
afterwards . 

M N M,N 

We first need to establish an upper bound on at ’ (v ’ ) . 

Let w^’^(x,t) = 0. Then 

( 5 . 3 . 8 ) * M ' N (« M ' N ) * || f 2M - 1 ' 2N ' 1 || 2 * || f 2M -‘ || 2 , , 

O.QxJ Qx{-1} 

* II F II 2 , + II f II 2 , ,, 

O.QxJ Qx{-1} 

using (5.3.2). Hence we can conclude that 

( 5 . 3 . 9 ) at M ' N (v M,N ) £ C. 

To estimate || L M * N v M ’ N - F 2 *" 1 * 2 *" 1 || we need to bound 

" -s.QxJ 
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(l M,N v M,N _ p2M-l , 2N-1 


II * II 


for <j> g H. 


s ,QxJ 

Consider the periodic IVP 


(5.3.10a) 


(5.3.10b) 


L ’ ip = <p for (x,t) e QxJ , 
ip = 0 for ( x , t ) e Qx { - 1 } . 


Then ip is a smooth function and using estimate (5.2.3) we have 

(5.3.11) || ' || i C || * || 

s.QxJ s,QxJ 

where C is a constant which depends only on the smoothness of 
the coeffecients of the modified IVP and hence of the original 


Let Q ’ be the projection operator that maps functions 


belonging to H f] H 


into V M ’ N defined as: 


1 , QxJ 


Q^ ,N w is the unique element of such that 


n ~M,N 

|| w - Q w 


l,QxJ 


s M,N V M,N 


|| w - s 


1 ,QxJ 


Then it is known that 


(5.3.12) 


_M,N ,, 
w - Q w || 


1 jQxJ 


s C N 1 S || w || 


s ,QxJ 


t 4 . ~ M,N _ n M,N 

Let ip = Q ip. Now 


( l M.N v M.N . p2M-l , 2N-1 ^ # J 


= ( l M ' N v M ’ N - f 2M - 1 ’ 2N ' 1 . l m - n v ) 

QxJ 

( l M,N v M,N _ jj2M-l , 2N-1 * l M,N ~ M,N } 


* ( l M.N v M,N . P 2M-1,2N-1 l M,N ( v _ ~ M.N,, 

QxJ 
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But by (5.3.7) 

( L M,N v M#N - f 2M " 1s2N_1 , l M,N y M ’ N ) 

+ ( V M S N _ ? 2M-1 t J M.N } = Q 

Ox { - 1 } 

Now since ip = 0 for (x,t) e Ox{-l) we may write 


OxJ 


< v M ' N - f 2 “-» , ; m ’ n > 




c v M ' N - 1™- 1 , ; - „ > 


flx{-l J 


Hence we can conclude that 

< l M - N v M,N - p2M-l ,2N-1 ^ , 

OxJ 

(K o i _ ( „M,N ? 2M-1 ~ M,N x 

(5.3.13) =(v -f , ip - ip ) 

Ox { 1 } 

+ ( l M ' N v M ' N - F 2M -‘- 2N - 1 ,L M ' N C V - ; M - N )> 
Now using (5.3.12) we can conclude that 


OxJ 


.M.N, MjNvii n . . 1 s || ii 

L (ip- ip )|| ^ C N II ip || 


0 ,OxJ 

And applying (5.3.11) we may write 


s ,QxJ 


(5.3.14) || L M ’%- i M,N ) || £ C N 1 S || <t> || 

0 , Ox J s , Ox J 


But 

(5.3.15) 


| ( l M,N y M,N _ f 2M-l,2N-l T M>N ( - tn _ ; M,N 


, L“*"0 p- ip “’")) | 

OxJ 


£ || L M ‘%- y> M,N ) || x || L MjN v M,N - F 2M 1,2N_1 || 


0 ,QxJ 

* C s N 1_S || * || 


0 


s ,OxJ 


using (5.3.9) and (5.3.14). 
Next , we estimate 


k V M,N - f 2 ”- 1 . ; M ' N - * > 


Ox { - 1 } 


,OxJ 
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From (5.3.9) we have that 


(5.3.16) 


|| v M ’ N - f 2 * 1 ' 1 || 


£ C. 


0 , Ox (-1 } 


li ~ M,N .. 

II v - V II 


0 ,Qx{-l } 


* O II J M,N - r II 


1 ,QxJ 


by the trace theorem, and so by (5.3.12) we obtain 


n 7 M,N .. 
II V ~ V || 


£ C N* 8 || V || 
0,fix{-l) s.HxJ 


Using estimate (5.3.11) once again we conclude that 


(5.3.17) 


ii 7 M»N ,, 

II V ~ V II 


* C N 1 8 || <p || 

0 , Qx {-1 } s.QxJ 


Hence applying (5.3.16) and(5.3.17) we get 


(5.3.18) 


M,N 7 2M-1 ~ M,N . 

I ( v ’ - f , ip -ip ) 


Qx{-n 


* C a N 1 - 3 || * || 


s.QxJ 


Combining (5.3.13) , (5.3.15) and (5.3.18) we obtain 


M,N v M,N _ f 2M-1,2N-1 | #) N l-s „ # , f 

Ox J S s , Qx J 

and this gives us the required estimate 


.3.19) || L M,N v M,N. ? 2M-1.2N-1 


-S.QxJ 


^ C N 
s 


Next, we estimate 


II v M ' N - ? 2M -‘ || 


-s ,Qx{-l } 


Consider the periodic IVP 


(5.3.20a) 

(5.3.20b) 


L M,w ip = 0 for (x , t ) e QxJ , 
ip = n for (x,t) e Ox { - 1 } 


Then ip is a smooth function and using estimate (5.2.3) we have 

(5.3.21) || V || s C || (i || ■ 

s ,QxJ s , Qx { - 1 } 

Let ip M,N = Q M,N ip. Now 
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( v M * N - f 2 ^ 1 i#l ) 


= c V M,N - m >n } 

Qx {-1 } Qx{-1 } 


+ ( v 


.M.N -2M-1 


, ip - ip 


M.N 


Qx { - 1 } 


But by (5.3.7) 

( L M,N v M,N - f 2M_1)2N_1 , l M jN i M,N ) 


ftxj 

+ ( y M,N _ f 2M-l § ; M.N 


) 


0. 


■ M.N 


Qx{-1) 


And since L ' ip = 0 for (x,t) e QxJ we may write 
( l M,N y M,N _ p2M-l , 2N-1 l M,N ~ M.N ) 


QxJ 


= C L M ' N v M ' N - f 2M-1.2N-l i l M,N ( - M.N . , 


ftxj 


Hence we can conclude that 


C v M * N - f 2 * 1 " 1 ) 


(5.3.22) 


Qx { - 1 } 

_ , M,N ? 2M-1 , ~ M.N , 

- ( V " f > * ~ V >Q X {-1} 

+ ( L M ’ N v M ’ N - f2M-1.2N-1 l M.N ( v| _ ~ M. N)) 

QxJ 


s: C N 
s 


1— s 


Thus we can show 
(5.3.23) || v M,N - f 2M_1 || 

-s ,Qx{-l } 

using (5.3.22) and the arguments employed earlier. 
We now need to estimate 


|| l M,N u _ p2M-l , 2N-1 


-s.QxJ 


We know that u satisfies u t - A u x - B u = F in the sense of 
distributions. Accordingly we may write 

(5.3.24) L M,N u - f 2M " 1,2N_1 = (L MjN u - Lu) - (f 2M_1 * 2N_1 - F) 

= -(a M “ 1,N-1 - A) u _(b M_1 )N_1 - B) u _(f 2M-1 * 2N_1 - F). 


x 


Now by (5.3.4) 
(5.3.25) || A - A 


M-l , N-l 


s.oo.QxJ 


* C s N _S || A || 

3s, <», QxJ 


and so 
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(5.3.26) || a M-l.N-1 || A C || A || 

s,a>,ftxj 3s.oo.ftxJ 

for M and N large enough. 

Let us show how to estimate the various terms in (5.3.24). It is 

known (Canuto at. a/. 1988) that the projection operator has the 

property that 

(5.3.27) || F - i? 2 * 1 " 1 ’ 2 **- 1 || * C N" S || F || 

“S.ftxJ O.ftxJ 

Next , we shall estimate || (B - §M~ 1 . N- 1 ) u |j 

-s.ftxJ 


For this we need a lemma. 


Lemma 5.3.1 Let A e H 


and v e H 


s ,00 ,ftxj 


. Then Av e H 


-s.ftxJ 


and 

(5.3.28) || Av || 

We have 


* C l| A || 


s " "" " x II v II 

s.QxJ s.oo.ftxJ -s.ftxJ 


Hence 


(Av, 4>) = (v, A 4>) , by definition, 

ftxj ftxj 


| ( Av , $) | | (v, A <J») | 

1 ftxj 1 _ 1 ftxj 1 


II ♦ II 


S.ftxJ 


II <0 II 

S.ftxJ 

I (v, A%) | || A% || 

ftxj s.ftxJ 


II * ♦ II II *11 , 

S.ftxJ s.ftxJ 


And this gives 

II Av || 


II A% || 


s.ftxJ 


£ sup 

(-s) ,QxJ ^ € g II ♦ II 


X V 


-s .ftxj 


-s .ftxj 


S.ftxJ 



Now it is easy to see that 


* . W 

ii * * ii 


S.QxJ 


sup 
0 € H 


II ♦ II 


* c s II A II 


s.oo.QxJ 1 


s , Qx J 

And this given us the required result. 

Thus we obtain 

j|(B - B M_1 ’ N_1 ) ulj <; 

-s.QxJ 

<= 8 ll» - b m_1 ' n_1 |I X || U II 

s.oo.OxJ -S.QxJ 


But 


and 


II u II * II u || 

“8, QxJ O.QxJ 


|| B - § M “ 1,N “ i || 


* C N“ S || B || 

8 ,oo ,OxJ 3 s,oo,QxJ 


Hence we obtain 


(5.3.29) || (B - B M_1 ,N-1 ) u || 


S ,oo,QxJ 


£ C N _S . 
s 


Next » we estimate || u j| 


Let 0 € H . Then 


■s ,QxJ 


Cu , 0) * - (u, 4> ) 

* Qx J x QxJ 

since both u and 0 are periodic in x. 

Hence 


j (u x . 0) | | (u, 0 x ) 


QxJ 1 


QxJ 


II ♦ II 


s.QxJ 


II * II 


s.QxJ 


But 


II ♦, II 


(s-l ) ,Qx J 


* II ♦ II 


s.QxJ 


And so we can conclude that 
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8 UP 


l (l V *> I 

1 * QxJ 


(u ,<j> x ) 


sup 


QxJ 


♦ e « II ♦ I! 


s ,Qx J 


♦ € H II * x II 


8 - 1 , fix J 


which gives us 

(5.3.30) || u x || * || u || 

~ H , Qx J -s+1 , QxJ 

But 

(5.3.31 ) || u || « || u || 

-s+1, QxJ O.QxJ 

And so by the lemma just proved we get 

(5.3.32) || (A - a M " 1,N “ 1 ) u || i C N“ 8 . 

s , co , QxJ 

Combining all these estimates we get the required result 

(5.3.33) |, l H.N „ . f »-1.2H-l „ s C N -8 . 

-s , QxJ 

Also we have (Canute et.al, 1988) 

(5.3.34) || u - f 2 **” 1 || i C N _s . 

— s , Qx ( — 1 } 

We can now prove our main theorem. 


Theorem 5.3.1 Let v M,N be the solution obtained by minimizing 
a H ' N <w M,N > as described in (5.3.5). Then for all s i 0 the 


estimate 


(5,3.35) || u - v M,N || 


■8 , Qx { 1 } 


+ || u - v H ' N || S C s N 

-s ,QxJ 


1-s 


holds . 

We have by (5.3.19) that 


j, l M,N v M,N _ f 2 M-l, 2 N-l || * c s N 

-s.QxJ 


1-s 


Moreover , by (5.3.33) we know that 
M,N f2M“ 1 , 2N-1 


II c 


11 4 C B N 

"8 ,QxJ 


1-s 


u 
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Using t hr triangle inequality we obtain 


(5.3. 3b) 1| L (u - v M > N ) || i C N 

-S.QxJ 8 

Finally , we have 


1-s 


(5.3.37) || u - v M,N || 

~B , OX ( - 1 } 

* II u - f 2M ~ ! || 

“8 , Ox { - 1 } 

£ C s N ! " s , 



V M,N 


II 

~s , Qx (-1 } 


using (5.3.23) and (5.3.34). 

Therefore using estimate (5.2.4) along with (5.3.36) and (5.3.37) 

we conclude that 


| u - v M,N || + || u 

s.OxU } 


M.N 


II 

-S.QxJ 


s C* N 

H 


l-s 


5.4 Hecover ing Pointwiae Values with Spectral Accuracy 

In this section we briefly describe how the local smoothing 
proposed by Gottlieb, Tadraor < 1 885) and Abarbanel , Gotti ieb(1985) 
and can be used to recover pointwiae values with spectral 
accuracy at any point in a neighbourhood of which the actual 
solution is smooth. If we wish to recover the values at t = 1 the 
local smoothing is particularly simple. Suppose we wish to obtain 
the value of the solution at the point (Xq,1). We assume that 
there exists a neighbourhood 

J = { t ; I x - I & 6 I 

in which the actual solution u(x , t) is smooth. Let p(x) be a C Q 
function with support in the set J and such that p is nonnegative 
everywhere and p( x Q ) * 1. Choose K = with 0 < 0 < 1 • and let 

It 

D (£) denote the Dirloblet kernel 
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I) K ( f, ) 


Y „ 8 i n( (2K+ 1 ) £/2) 

J*-K «ln(f,/2) 


5 * 2mn , 


• 21 * 1 


. § =2mjr. 


Than to oh! am the regularized 
def i nc 


version of v M * N at ( X()) l) we 


,2k 


(5.4.1 ) K v * <x 0 .l> * **-*- J D K (x 0 - x) p( x ) v MlN (x,l) dx. 

2* n 


It has been proved in (Canuto et.al. 1988) that if 

H u - v M,N j) s C M“ s+1 . then 

-a.OxUJ 8 

(5.4.2) lu(* 0 , 1 ) - H v H,N (x 0 ,l)l s Cj (1 + log M ) M~ s+1 

* C 2 lT 8 *^ 8 , 

where t he eonutant n V and t* 2 depend upon the Sobolev norms of p 
and u mrr the interval J . A balance of the errors is aoheived by 

putting § « 1/2 , in which case we obtain 

1 U ( Xq , 1 ) - R v M,N (x 0 ,l)l « 0 ( M" s/2+1 ) , 

which proven that u(x () .l> can be approximated with spectral 

accuracy start mg from the knowledge of the Galerkin-Col location 

M N 

approximation v * . 

Suppose now that we wish to recover the value of the 
solution at an interior point (x 0 ,t Q ). We assume that u(x,t) is 
smooth in the set 0 , where 

0 * { <x,t) : lx — Xq I i 5 i It - t Q l s £ ). 

Let p(x) be a C® function with support in the set 

J * { x : lx - x^f < $ ) * 

which is nonnegative everywhere and such that p(* 0 ) 
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Sim? . nt J < 


* '» ».th .upper 

' " - > > ) . 


1 tn set 


C ) 5 t » Choose 


in nonnaga t j v*» ever ywhrr r *54 a _ , , 

" yU ‘ g rj(t 0 ) 3 *• Choose 

K * *t #i ami I * hi 1 « j t h o < « , . . , v 

' U ' ° "■> -•».». the 

” ! “ f k ' rn - 1 * m ’ ''"V <ha Legendre kerne, 

i 1 

' " • Z a '* * >'*> S<«. 


Then to obtain thr ra*u lari red v « | Mc ,„ 


define 


of v"'" at fv ♦ n 

W we 


* V " 0 '" ‘ ■i* IJj ’■'"•'o’ *><*> r,U) v M ' M {x,t>d*dt.' 

Once •or» t it mr% !»*» shown that k v^*^r* , . 

lx 0 ,l 0 ? approximates 

u(. 0 ,t 0 ) ».,h .por.re, .rrur.cy .ml opt,.., bp , npoe of the 

error* t « obt * t ni»d by choosing $ • , *1/3, 


1.5 Nu**r i c* J Result* 


I" ,h ‘" "" c,,an »• -»«°».»r.l. the erneecy of the eethod 

proposed in this paper 

Eianpl# 5 . 5 ,t 

Consider t ha init t*i value probles with periodic boundary . 
u i “ «<».n v M - b<*,t) u • F(x 1 1 ) 


0***2*. -1 I ti 1 


ith initial condition 


U(*,-| ) • g(«) 

•ml let if' » « h*s n diaeont inuty to tig derivative 
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and with same a(x,t> and b(x,t) . 

Results of smoothing of the spectral approximation of U(x,t) , 
with M = 128 and N = 17 are displayed below 
Case I 


X P = £ (p+1/2) 

v equals 

IUCx^.1) - V m ’ n (x pJ l)| 

|U(x y ,l) - R V m,n (x yJ l)j 

4 

1.47 (-3) 

2.69 (-8) 

5 

1.88 (-3) 

2.28 (-8) 

6 

2.32 (-3) 

2.67 (-8) 


Table 5.5.1 

Case II 


x p = \ 

v equals 

|U(x p ,l) - V m,n (x y ,l)| 

1 UCx , 1 ) - R V m,n (x y ,l)l 

4 

1.12 (-3) 

2.88 (-8) 

5 

1.11 (-3) 

4.19 (-8) 

6 

1.28 (-3) 

4.48 (-8) 


Tabl e 5.5.2 

where the number in bracket denote the power to base ten. 
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